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Hie  autonomous  firefighter  manipulation  problem  has  two  main  areas:  design 
of  a  manipulator  mechanism  and  controller  design  for  that  mechanism.  In 
addressing  these  problem  areas  the  following  assumptions  were  made: 

1 .  The  manipulator  must  manipulate  the  USAF  tool  shown  in  Figure  1 . 

2.  The  manipulator  must  be  as  autonomous  as  possible. 

3.  Due  to  (2),  the  raaniplator  must  be  as   power  efficient  as  possible  in 
order  to  conserve  on- board  power  supplies. 

4.  The  tool  positioning  requirement  would  be  approximately  +  0.10  in  and 
+  5  deg  at  the  drill  tip. 

The  mechanism  design  must  first  be  concerned  with  the  expected  loads 
associated  with  the  drilling  process.  Since  these  are,  in  general,  widely 
variable  depending  on  the  type  of  drill,  type  of  material  to  be  drilled, 
incidence  angle,  reinforcements,  and  many  other  parameters,  we  were  not 
surprized  to  find  very  litle  published  material  on  the  drilling  process  in  the 
open  literature.  To  remedy  this  we  undertook  several  experimental  projects  to 
measure  and  characterize  the  particular  process  of  drilling  aircraft  skins. 
Given  the  measured  drilling  load  characteristics,  the  manipulator  mechanism 
and  actuators  could  be  designed  or  selected. 

The  manipulator  controller  is  concerned  with  guiding  the  manipulator  from 
point-to-point  in  a  power  efficient  manner,  and  maintaining  a  tool  orientation 
while  a  hole  is  drilled.  The  following  discussion  is  about  the  first  part  of 
the  controller  problem,  the  point-to-point  control  and  the  tasks  associated 
with  solving  this  problem.  Clearly,  the  manipulation  requirements  outlined  in 
assumption  4  above  could  not  be  considered  high  performance  specifications. 
Instead,  the  controller  of  interest  was  one  aimed  at  effective  coordination  of 
movement  and  was  therefore  called  a  low  performance  controller. 


The  remainder  of  this  report  is  divided  into  three  sections:  the  next 
section  sunmarizes  the  major  results  acheived  to  data;  the  following  section 
describes  the  most  recent  mainpulator  controller  development  project;  the  last 
section  discusses  final  drilling  load  measurement  studies. 
II.  Results 

A.  Facilities  Developed 

Perhaps  the  most  important  facility  that  has  been  developed  for  this 
project  is  the  low  Performance  Control  (IPC)  lab.  This  test  bed  facility 
consists  of  a  computer  based  controller  system  designed  to  be  used  to 
investigate  the  merits  of  various  candidate  controller  and  sensor 
philosophies  for  low  performance  plants.  The  next  section  of  this  report 
describes  the  IPC  lab  in  more  detail  when  coupled  to  a  manipulator  so  it 
will  only  be  mentioned  briefly  here.  The  lab  presently  consists  of  an  IBM 
PC- AT  connected  to  a  rigid,  revolute  robot  (Neptune  by  Alan  Hayes  Corp.). 
The  Computer/manipulator  interaction  is  accomplished  via  LABPAC  hardware 
and  software.  In  addition,  hardware  for  optical  encoding  of  joint  data, 
and  an  improved  data  taking  and  display  system,  including  direct  memory 
access,  will  be  purchased  and  installed  in  the  next  few  months.  This 
facility  has  been,  and  will  be,  a  valuable  and  necessary  tool  for  research 
into  computer  based  control  of  low  performance  machinery  at  the  Naval 
Postgraduate  School. 

The  IPC  testbed  is  shown  in  Figure  2.  The  Controller  (IBM  PC-AT)  is 
at  the  left  in  the  figure,  the  Data  System  (IBM  PC-XT)  is  at  the  right, 
and  the  external  patch  rack  is  in  the  center.  The  Data  system  can  be  run 
independently  of  the  controller  to  generate  plant  modeling  data,  or  the 
Data  System  and  the  controller  can  be  run  simultaneously  through  a  XENIX 


Figure  2.  Computers  For  the  Low  Performance  Control  Lab 
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operating  system.  Maximum  bandwidth  is  about  100  Hz  for  single  input 
single  output  operation. 

Ihe  Neptune  manipulator  is  shown  in  Figure  3  mounted  on  its  support 
stand.  The  hydraulic  pump  is  to  the  right  of  the  stand,  the  hydraulic 
accumulator  is  beneath  the  stand,  and  the  bottom  shelf  holds  the  valves 
and  a  hydraulic  manifold.  The  manipulator  has  six  revolute  degrees  of 
freedom,  operates  on  a  water  based  fluid  at  120  psig,  and  can  be 
configured  for  solenoid  ar  servovalue  operation. 

Besides  the  IPC  lab,  two  more  or  less  temporary  experimental 
facilities  have  developed  out  of  the  NPS  Robot  Firefighter  Project,  these 
are  the  drilling  apparatus  and  the  flexible  link  facility.  Key  features 
of  the  drilling  apparatus  were  shown  in  Figure  1 .  The  LEAF  drill  was 
mounted  on  an  "H"  mounting  and  afixed  to  a  movable  milling  bed  platform. 
The  vertical  support  in  the  "H"  mounting  was  a  steel  cylinder  of  1  inch 
diameter  which  was  instrumented  for  strain  in  order  to  estimate  the  loads 
incurred  by  the  drillnig  process.  The  use  and  analysis  of  the  drilling 
apparatus  is  described  in  detail  in  Section  IV  of  this  report. 

The  flexible  link  facility  is  shown  in  Figure  4.  The  facility 
consists  of  a  hydraulic  pressure  source  (2000  psig,  shown  at  left) ,  a 
support  platform  and  control  box  shown  at  the  center,  a  one  degree  of 
freedom  hydraulic  motor  actuator  at  right,  and  a  flexible  box  beam  shown 
vertically  upwards.  The  beam  was  designed  to  flex  in  one  plane  only  in 
order  to  study  distributed  dynamic  structural  flexibility.  This  facility 
was  designed  and  implemented  based  upon  our  desire  to  use  the  firefighter 
mechanism  to  accomplish  a  range  of  tasks- expecially  the  heavy  lifting 
associated  with  manipulating  bombs,  rockets,  or  large  aircraft  parts. 
Consequently,  if  the  mechanism  is  not  to  become  ridiculously  bulky  and 
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heavy,  it  can  be  expected  to  be  flexible;  that  is,  it  should  be  designed 
to  be  flexible.  Such  a  flexible  mechanism  poses  interesting  problems  in 
terms  of  modelling  and  control.  Much  of  what  we  have  learned  in  the 
flexible  link  studies  is  being  incorporated  into  space-based  manipulator 
analysis  at  the  time  of  this  writing. 
B.  Major  Conclusions 

The  conclusions  which  follow  are  referenced  to  the  theses  in  which 
they  were  first  mentioned.  The  theses  are  listed  in  the  next  portion  of 
this  report. 

1 .  Manipulator  Design  for  the  USAF  Drill 

(a)  Drill  walking  at  the  inception  of  drilling  can  lead  to  large 
lateral  forces  and  corresponding  twisting  moments  being  exerted  on 
the  manipulator  support  structure  (Lawrence) .  This  can  be  largely 
avoided  by  pilot  hole  drilling  or  drill  tip  redesign. 

(b)  Manipulator  load  forces  vary  significnatly  as  the  drill  tip 
wears  (Lawrence) .  Design  of  the  manipulator  must  account  for  the 
worst  case  to  be  expected.  At  some  point  the  drill  becomes  so 
worn  that  it  is  unable  to  penetrate  the  skin  (Yobs) .  A  means  for 
sensing  the  onset  of  this  condition  should  be  provided  in  the 
system  design  and  utilization. 

(c)  The  drill  tip  loads  have  been  characterized  for  both  the  worn 
drill  (Yobs)  and  the  new  drill  (Lawrence) .  Dynamic  and  peak  loads 
have  been  measured.  This  data  will  allow  the  subsequent  design  of 
mechanism  and  acutator  components  for  a  drill  manipulator. 

(d)  Flexibility  effects  are  difficult  to  predict  and  measure 
(Petroka) .  This  conclusion  suggests,  that  two  manipulators  may  be 
a  good  idea:  one  being  ridid,  low  powered,  and  low  performance  to 


manipulate  lightweight  tools  such  as  the  USAF  drill;  and  one  being 
strong,  high  powered,  and  perhaps  high  performance  to  manipulate 
heavy  objects.  The  strong  manipulator  may  be  flexible  and  would, 
consequently,  require  a  very  sophisticated  control  system,  perhaps 
even  a  vis ion- based  controller. 

2.  Rigid  Manipulator  Modeling 

(a)  Simulation  of  robot  motions  should  be  done  in  global 
coordinates  to  avoid  singularity  problems  associated  with  local, 
or  joint  coordinate  formulations  (McCarthy) .  Such  global 
simulations  are  easier  to  formulate  and  understand  than  local, 
recursive  forms  (Mohammed) . 

(b)  For  motions  less  than  30  degrees,  the  Neptune  robot  appears  to 
be  a  very  nearly  linear  machine  (Harris ,  Lewis) .  This  result  has 
strong  implications  in  the  design  of  controllers  for  low  speed 
manipulators  such  as  the  Neptune. 

3.  Rigid  Manipulator  Control 

(a)  The  USAF  drill  should  not  be  used  at  incidence  angles  greater 
than  about  12  degrees  from  the  local  normal,  or  at  locations  which 
are  supported  by  stiffners  beneath  the  skin  (Yobs) . 

(b)  Optimal  Control  (LQR)  has  been  demonstrated  as  an  appropriate 
means  of  achieving  effective  coordinated  control  of  rigid  link 
motions  (Harris) . 

C.  Thesis  Advised 

1.  McGalliard,  G.  R. ,  "A  General  Simulation  Program  for  Robot 
Manipulator  Arm  Dynamics,"  NPS  MSME  Thesis,  Sept.  1984. 

2.  MaCarthy,  W.  F.,  "Simulation  of  High  Speed  Motion  of  Rigid, 
Revolute  Mechanisms,"  NPS  MSME  Thesis,  Dec.  1985. 
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3.  Lewis,  D.  R. ,  "Modelling  of  a  Low  Performance  Rigid  Revolute  Robot 
Arm,"  NPS  MSME  Thesis,  Dec.  1985. 

4.  Yobs,  R.  L. ,  "Manipulator  Load  Forces  for  a  Robotic  Firefighter," 
NPS  MSME  Thesis,  Mar.  1986. 

5.  Petroka,  R.  P.,  "Computer  Simulation  and  Experimental  Validation 
of  a  Dynamic  Model  (Equivalent  Rigid  Link  System)  on  a  Single- Link 
Flexible  Manipulator,"  NPS  MSME  Thesis,  June  1986. 

6.  Burrill,  L.  D. ,  "A  Feasibility  Investigation  for  Optimal  Robotic 
Control,"  NPS  MSME  Thesis,  Sept.  1986. 

7.  Mohammed,  K. ,  "Non- Singular  Modelling  of  Rigid  Manipulators,"  NPS 
MSME  Thesis,  Dec.  1986. 

8.  Harris,  J.  P.,  "Investigation  and  Development  of  a  Micro- computer 
Based  Robotic  Controller,"  NPS  MSME  Thesis,  June  1987. 

D.  Theses  in  Progress 

1 .  Sanders ,  D. ,  "Optimal  Control  as  an  Autonomous  Vehicle  Path 
Planner,"  NPS  MSME  Thesis,  Dec.  1987  (expected). 

2.  Fancher,  C,  "Microcomputer  Based  Optimal  Control  for  a  Low 
Performance  Manipulator,"  NPS  MSME  Thesis,  Dec.  1987  (expected). 

3.  Altinok,  S.,  "Non- Singular  Simulation  of  Rigid  Manipulator 
Motions,"  NPS  MSME  Thesis,  Dec.  1987  (expected). 
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T,Ti,T2  — ■  Joint  torques 

u  --  Input  matrix 

Ui,LT2  —  System  inputs 

v,vi,V2  —  Servovalve  input  voltages 

x  --  State  matrix 

Pl,p2  -  Joint  friction  coefficients 

*  --  Discrete  plant  matrix 

r,f  0,1*1  —  Discrete  plant  matrices 

A.  --  Scalar  factor 

0.  -  Servovalve  electrical  resistance 

§>§1>§2  -  Joint  positions 

0,01,02  —  Joint  velocities 

x  --  Time.delay 
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situations.  Deep  submergence  undersea  exploration  vehicles,  unexploded 
munitions  disposal  units  and  carrier  flight  deck  firefighting  have  stood  out  as 
examples  of  jobs  that  could  be  taken  over  by  robots  with  the  goal  of  reducing 
the  risk  to  human  life  [Refs.  1,2,3].  In  the  area  of  shipboard  fires,  it  was 
realized  that  a  motorized  firefighting  unit  cannot  replace  a  complete 
firefighting  team,  but  a  quick  reaction  vehicle  could  provide  intial  fire 
suppression  and  quite  possibly  keep  the  fire  from  spreading  out  of  control  [Ref. 
11.  Additionally,  the  vehicle  could  be  sent  into  areas  where  the  fear  of  either 
ammunition  or  fuel  explosion  would  prohibit  a  firefighting  team  from 
entering. 

The  design  and  operation  of  a  control  system  for  such  a  robotic  device  is 
one  area  that  has  recently  come  under  renewed  scrutiny  [Ref.  2].  The  majority 
of  all  robots  available  today  use  a  somewhat  large  and  immobile  computer  as 
the  master  "brain".  These  large  computers  are  needed  to  handle  voluminous 
calculations  that  are  demanded  by  the  complex  control  algorithms  for 
multiple  robot  manipulator  links.  In  light  of  their  sheer  size  however,  these 
systems  are  a  poor  choice  for  mobile  operations.  A  more  suitable  choice  is  a 
micro-computer  based  controller.  Due  primarily  to  increases  in  speed  and 
computational  power,  "personal  computers"  (PC's)  and  the  micro-processors 
that  drive  them,  have  been  thrust  into  the  limelight  as  the  bridge  to  the  self- 
contained  robot  operations. 

The  first  goal  of  controller  development  is  the  selection  of  the  algorithm 
that  will  provide  adequate  system  performance.  Because  cost  is  a  major  factor 
in  most  operations,  an  ideally  suited  controller  should  minimize  the  power 
required  to  conduct  an  operation  while  balancing      this      against    the 
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I.  INTRODUCTION 

Robotics,  only  recently  considered  science  fiction,  has  now  become  a 
permanent  part  of  technology.  Although  the  term  "robot"  is  often 
misinterpreted,  robotic  type  devices  can  now  be  found  in  all  aspects  of  society, 
from  personal  toys,  to  an  integral  part  of  many  major  manufacturing 
industries.  For  the  purpose  of  this  paper,  a  robot  is  defined  as  a  self  contained 
electro-mechanical  system  consisting  of:  a  digital  controller,  a  power  source, 
and  a  manipulator  arm. 

The  advantages  of  using  machines  to  replace  man  -have  long  been 
appreciated.  In  industry,  economics  has  driven  the  incorporation  of  this  new 
technology  to  accomodate  and  sometimes  replace  older  methods.  The 
resulting  automation  has  markedly  boosted  output  while  decreasing  the  per 
unit  manpower  costs.  This  has  in  turn  caused  an  overall  increase  in 
production  efficiency  making  many  items  more  affordable  to  the  consumer. 
Another  advantage  of  automation  has  been  the  attainment  of  more  consistent 
results.  Painted  items  thus  become  more  uniform,  machined  items  become 
more  exact,  and  production  quality  becomes  more  standard  from  item  to  item. 
All  of  these  facts  not  withstanding,  however,  the  most  important  reason  for 
the  use  of  mechanical  equipment  is  to  lessen  the  endangerment  to  man  while 
operating  in  hazardous  environments.  The  more  autonomous  the  robotic 
devices,  the  more  removed  and  hence  safer  the  person  controlling  the 
mechanism. 

Recently,  the  U.S.  Navy  has  developed  an  interest  in  the  use  of  robots 
with  an  aim  toward  replacing  humans  in  potentially  life  threatening 


r. 
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positional  error  generated  by  the  control  system.  To  accomplish  this,  a 
modern  controller  should  be  investigated  for  implementation.  Such  modern 
controllers  are  capable  of  coordinating  multiple  joint  inputs  to  control  the 
complicated  manipulator  output  configuration.  The  proof  of  these  controllers 
lies  in  their  ability  to  control  a  real  piece  of  hardware. 

This  thesis  describes  the  implementation  of  a  modern  controller  on  the 
Naval  Postgraduate  School  manipulator  control  test  bed.  The  test  bed  consists 
of:  a  NEPTUNE  multiple  rigid  link  robot  system  (Fig.  1),  an  IBM  PC-AT 
micro-computer,  and  two  Atchley  jet  pipe  electro-hydraulic  servovalves.  The 
NEPTUNE  system  was  modified  so  that  it  could  be  run  as  originally 
purchased  (solenoid  operation),  or  through  a  PC  based  servo  controller  in 
order  to  demonstrate  the  comparative  performance  associated  with  each  type 
of  controller. 
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Figure  1 
NEPTUNE  Manipulator 
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II.    BACKGROUND/PROBLEM  DESCRIPTION 

The  primary  goal  of  this  work  was  to  initiate  an  NPS  test  bed  system  for 
Navy  development  of  robotic  controllers.  The  fact  that  operational  robots  will 
be  needed  to  operate  in  a  hostile  environment  suggests  that  water  based 
actuators  be  used  as  the  prime  motive  force.  The  NEPTUNE  robot  arm  and 
control  system  fit  this  requirement  and  were  previously  selected  [Ref.  4]as  the 
system  nucleus.  Lewis  [Ref.4]  described  the  operation  of  the  basic 
computer/solenoid  controlled  assembly  which  was  modified  in  the  present 
investigation. 

In  the  original  NEPTUNE  actuation  scheme,  a  pseudo-proportional 
control  was  obtained  through  the  use  of  several  fluid  restrictors  (Fig.  2).  By 
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Figure  2 
Solenoid  hydraulic  system 
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ganging  more  than  one  of  these  restrictors  together  on  a  joint,  the  motion  and 
flow  rate  history  as  seen  in  Figure  3  was  acheived.  However,  in  order  to 
investigate  advanced  servo  controls,  an  infinitely  variable  flow  rate  is 
required  as  in  Figure  4.  In  order  to  modify  the  NEPTUNE  to  accept  a  servo 
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Figure  3 
Solenoid  operation,  (a)  motion,  (b)  flowrate  history. 
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Figure  4 
Servo  valve  operation,  (a)  motion,  (b)  flowrate  history 
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control  signal,  Burrill  [Ref.  5]  conducted  a  feasibility  study  to  determine  the 
plant  requirements  and  select  the  best  suited  servo  valve.  "Off  the  shelf' 
valves  were  desired  to  both  minimize  the  initial  procurement  costs,  and  also  to 
standardize  the  system  components  with  those  generally  available  to 
industry.  The  Atchley  model  218  jet  pipe  servo  valve  was  selected,  and  two 
valves  were  placed  in  service.  Because  only  two  dimensional  control  was  to  be 
demonstrated,  axes  1  and  2  were  chosen  (Fig.  5). 


AXIS  S  270 


AXIS  0 


Figure  5 
NEPTUNE  axes  of  rotation 


An  initial  hardware  modification  was  required  to  improve  hydraulic 
system  operation.  To  provide  a  reservoir  of  fluid  energy,  the  NEPTUNE 
hydraulic  power  pack  comes  with  a  42  cubic  inch  air  bladder  accumulator.  In 
order  to  reduce  entrained  air  in  the  water-based  operating  fluid,  an  air 
bladder  in  the  accumulator  is  filled  to  approximately  55  psi,  or  one  half  of 
system  operating  pressure.   The  fluid  pump   then  fills  the  accumulator  and 
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shuts  off  at  110  psi.  As  fluid  is  drawn  down  in  an  arm  move,  an  unloading 
valve  monitors  the  pressure  and  automatically  turns  the  pump  on  at 
approximately  100  psi.  With  two  joints  being  operated  simultaneously,  pump 
pressure  waves  can  be  observed  in  the  actuator  cylinder  data.  This  effect  was 
difficult  to  predict  and  evaluate. 

To  alleviate  this  problem,  joints  one  and  two  cylinder  dimensional 
measurements  were  taken,  and  the  maximum  amount  of  fluid  needed  to  move 
both  actuators  was  calculated.  It  was  determined  that  approximately  36  cubic 
inches  of  fluid  was  needed  to  move  both  actuators  the  length  of  their  travel.  A 
580  cubic  inch  accumulator  was  fitted  such  that  less  than  10  percent  of  the 
total  volume  would  be  used  on  any  arbitrary  move.  The  resulting  cylinder 
pressure  data  was  much  more  well  behaved. 

At  this  point  some  clarifications  should  be  made  as  to  the  distinction 
between  the  plant  and  controller  as  described  by  this  research.  Figure  6 
shows  the  general  multi-port  layout  with  the  term  system  defined  as  the 
combination  of  the  computer,  both  servovalves,  and  their  associated  links. 
The  term  plant  refers  to  the  arm  and  servovalve  combination.  Although  not  a 
radical  departure  from  conventional  thinking,  the  inclusion  of  the  servovalve 
in  the  plant  block  instead  of  the  controller  block  could  be  considered  unusual. 
There  were  two  reasons  for  taking  this  approach.  First,  the  plant  is  normally 
considered  as  that  part  of  the  system  that  is,  for  one  reason  or  another,  already 
selected,  unchangeable  and  must  be  controlled.  Secondly,  as  will  be  discussed 
later,  the  dynamic  characteristics  of  the  links  when  combined  with  the 
hydraulic  characteristics  of  the  valves  allows  a  simple  yet  effective  modeling 
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of  the  plant.  Because  the  speed  and  computational  efficiency  of  the  IBM  PC 
AT  are  significantly  less  than  that  of  a  larger  mainframe  computer,  any 
simplification  of  the  model  may  allow  for  a  shorter  controller  sample  rate.The 
controller  is  defined  as  the  combination  the  sensor  acquisition/conversion 
elements,  the  IBM  PC  AT  computer,  and  the  control  program.  With  the 
exception  of  Digital-to-Analog  Converters  (DAC),  Analog-to-Digital 
Converters  (ADC),  and  the  computer  itself,  all  control  and  regulation  is 
software  based.  The  DAC  and  ADC  are  contained  on  two  expansion  port 
boards  internal  to  the  AT.  Computer  interfacing  was  done  with  the  system 
known  as  LABPAC,  which  has  6  12  bit  DAC's  and  16  12  bit  ADC's.  With  12 
bits  of  data  availabe  for  conversion,  the  maximum  available  resolution  is  .024 
%.  In  addition  to  the  converters,  LABPAC  has  an  onboard  clock  timer  for 
sample  rate  regulation.  Once  setup,  LABPAC  automatically  does  both  data 
acquisition  and  control  at  set  intervals. 

The  system  model  is  based  on  the  Multiple  Input  Multiple  Output 
(MIMO)  elements  of:  the  computer  controller,  servovalves,  hydraulic 
actuators,  and  the  arm  itself.  To  create  motion  in  the  system,  a  voltage 
potential  (v)  is  applied  to  the  servovalve  electrical  input,  which  in  turn  sets  up 
a  current  flow  through  the  valve  actuating  coil.  This  causes  an  opening  of  the 
valve  orifice.  The  actuating  voltage  in  combination  with  the  actuator 
feedback  of  fluid  flowrate  (q)  produces  the  valve  pressure  outputs  Pi  and  P2. 
In  turn  Pi  and  P2  are  used  as  inputs  to  the  actuator  along  with  the  arm 
velocity  (0)  in  order  to  generate  the  torque  (T)  which  creates  motion  in  the 
arm. 
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In  order  to  validate  the  robot  test  bed  system  for  future  operation,  an 
optimal  state  space  controller  (LQR)  was  developed  and  implemented.  As 
discussed  by  Ogata  [Ref.  6]  and  Owens  [Ref.  7],  optimal  control  theory 
requires  that  a  merit  function  J  be  selected.  J  can  be  considered  the 
performance  index  that  is  to  be  minimized.  To  meet  with  established 
conventions,  the  merit  function  usually  uses  the  two  quadratic  forms  that 
follow. 

J^  JJ(xTQx)  dt  (Eqn.  I) 

J2~  %  (uTRu)  dt  (E(in-  2) 

The  matrices  x  and  u  are  the  state  space  matrices  that  correspond  to  state 
error  and  system-inputs  respectively.  Q  is  a  weighting  matrix  that  determines 
the  relative  importance  of  each  term  in  the  state  matrix.  R  is  determined  by 
the  system  performance  factor  to  be  minimized,  and  cost  weighting  given  each 
input.  J  can  take  the  form  of  the  arithmetic  combination  of  Equations  L  and  2 
in  the  form  : 

•/  =  •/.  -!-./.,  (Eqn.  3) 

J=  Jg{\(xrQx)  +  (urRu)}cfc  (Eqn.  4) 

This  is  helpful  because  for  most  real  life  systems  the  minimization  of  the 
state  error  alone  must  be  tempered  by  at  least  two  considerations.  First,  all 
operations  "cost"  something  whether  it  is  energy,  time,  or  some  other 
quantity.  Secondly,  since  this  form  of  optimal  control  produces  the  feedback 
law: 
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u  =  -Kx  (Eqn.5) 

a  large  state  error  will  consequently  produce  a  correspondingly  large 
correcting  input.  Too  large  a  feedback,  however,  can  saturate  most  real 
systems  thereby  giving  them  a  high  degree  of  non-linearly. 

To  narrow  the  scope  of  the  problem  and  provide  a  realistic  goal,  the 
performance  indices  chosen  were  related  to  state  error  versus  the  square  of  the 
hydraulic  power  used  by  the  system.  For  the  purpose  of  this  discussion,  the 
definition  of  hydraulic  power  will  be  that  used  by  Merritt  [Ref.  8]  where: 

Power  =P    q  (Eqn.6) 

P    =  P       ,    =  constant 
s  supply 

Initial  dynamic  modeling  of  the  two  link  arm  mechanism  was  done  in 
global  vice  local  coordinate  systems.  This  decision  was  made  in  anticipation  of 
future  problems  with  solving  for  joint  accelerations  using  the  Newton-Euler 
forward  dynamic  equations.  If  local  coordinates  are  used  in  conjunction  with 
the  Denavit-Hartenburg  [Ref.  9]  convention,  the  potential  exists  for  a 
singular  solution  when  two  adjoining  links  are  colinear  with  the  relative 
angle  being  either  0  or  180.  degrees  It  was  realized  that  the  joint  solutions  in 
local  coordinate  systems  are  more  computationally  efficient;  However,  the 
singularity  condition  described  above  precludes  general  implementation  of  a 
model  reference  position  control  algorithm  that  does  not  have  error  handling 
procedures. 
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Subsequent  modeling  of  the  NEPTUNE  system  revealed  that  a  much  less 
complex  plant  model  was  needed  to  satisfactorily  describe  the  motion.  This 
realization,  when  coupled  with  the  fact  that  all  sensors  and  actuators 
functioned  in  local  coordinates,  allowed  final  implementation  of  the  controller 
in  joint  relative  coordinates. 

There  have  been  numerous  research  efforts  aimed  at  various  types  of 
micro-computer  robotic  control  [Refs.  10-14].  Optimal  path  control  is  dealt 
with  by  Snyder  [Ref.  15].  The  closest  research  done  in  the  area  of  optimizing 
power  appears  in  Reference  16,  where  performance  indices  similar  to  the  one 
mentioned  above  are  used,  with  the  focus  on  weighted  near  minimum  time- 
fuel  control. 
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m.   ANALYTICAL  SYSTEM  CHARACTERIZATION 

A.    RIGID  TWO  LINK  SYSTEM  MODEL 

The  Newton-Euler  formulation  for  the  robot  arm  was  used  in  developing 

the  classical  dynamic  model.  Previous  work  by  Lewis  [Ref.  4]  on  the 

NEPTUNE  arm  showed  that  with  the  original  solenoid  valves  installed, 

reduced  dynamic  Equations  7  and  8  could  be  used  (all  angles  are  relative). 

Tl  =  Pl©2  +  D^sin(02~0l)"HG^sin(0l)"  Ms  (Eqn.  7) 

T2  =  P202  +  Dgsin(e2  -  0t)  (Eqn.  8) 

D  =m2d2+  rn3[l2  +  dj 

G  =  I  (m.,  *•  m  )  -i-  m   d 

\f    -  return  spring  moment 

Because  the  modified  arm  system  contained  two  servo  valves  of  a  higher 
performance  and  flow  rate,  the  full  Newton-Euler  two  degree  of  freedom 
equations  of  motion  were  used  [Ref.  17].  The  global  formulation  is  shown  in 
Figure  7. 

The  torque  required  for  each  joint  was  solved  for  as  follows  : 


Tl  =  (Ji  +  mld22+m2l2l)Ql  +  m2lld2cos(Ql-Q2)Q, 


2 


+  m2lld2sin(@l-e2)Q22   +(m.2g  l{  +  m{g d{)cosQl  +  T2  (Eqn.  13) 
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r2  =  (J2+m24)e2+m2/1cf2coS(e2-91)e1 


L2 


+  m2/l(f2sm(02-01)ei    +  m2gd2cos&2 


(Eqn.  14) 
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Figure  7 
NEPTUNE  geometry 


B.     SERVOVALVE  MODEL 

The  servo  valve  was  modeled  using  the  classical  equation  for  flowrate: 


a  =  ki\  P    -P, 

1  3  1 


(Eqn.  15) 


v 


Since  all  inputs  to  the  actuating  coil  were  DC,  the  valve  impedance  was 
modeled  as  a  simple  coil  resistance  (Q).  This  allowed  a  direct  conversion  from 
the  model,  in  current  (i)  to  actuating  voltage  (v)  (Eqn.  16).  With  the  servo 
valve  installed  in  the  system,the  flowrate  relationships  in  Figure  8  could  be 
used  determine  either  the  flow  or  pressure  drops  at  any  point  in  the  system,  if 
other  values  are  known. 
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q  =  kvy/P    -P, 

^  3  I 


(Eqn.  16) 


Pi 


P2 


Pr 


Figure  8 
Pressure/Flow  Schematic 


If  the  valve  is  assumed  to  be  symmetrical  (both  ports  and  spool  areas  are 
identical),  then  the  valve  constant  k  will  be  the  same  for  both  supply  and 
return.  Once  k  is  determined  experimentally  ,  Equation  16  can  be  rearranged 
as  follows : 


1  *  V  k  V  J 


(Eqn.  17 


Refering  to  Figure  8  and  Figure  2,  it  can  be  shown  that  0  and  the  flowrate  q 
are  linearly  related  by: 


q  =  c0 


c  =  .1131 


in 


deg 


{Eqn.  18) 
(Eqn.  19) 
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because  of  the  combination  piston,  and  rack  and  pinion  gear  assembly.  With 
P3  =  constant,  Equation  17  can  be  further  reduced  to: 


1         ^ch/ 


{Eqn.  20) 


or 


P    =fcniv,Q)  (Eqn.  21) 

C.    VALVE  TESTING 

Pressure/flowrate  tests  were  conducted  on  both  servovalves  in  order  to 

determine  the  valve  constant  k  (Table  I).  As  can  clearly  be  seen,  although  k 

did  follow  a  trend  for  each  run,  k  was  generally  not  a  constant-either  between 

the  ports  of  the  same  valve,  or  with  different  input  voltages.  One  explanation 

for  this  is  that  the  normal  operating  pressure  range  for  this  valve  is  between 

200-3000  psi.  The  pump  and  accumulator  currently  operate  at  between  100- 

110  psi.   It  is  possible  that  the  present  configuration  of  the  system  forces  the 

valve  to  operate  in  a  non-linear  region  below  the  minimum  design  point.  This 

notion  is  bolstered  by  the  fact  that  although  the  rated  maximum  input  voltage 

is  8.0  VDC,  the  observed  saturation  point  is  approximately  1.0  VDC  for  joint 

1,  and  0.6  VDC  for  joint  2.    At  inputs  greater  than  these  no  corresponding 

increase  in  flowrate  occurs. 

Pressure  plots  of  the  driving  cylinder  (Pj),  and  the  opposing  cylinder  (P2) 

showed  expected  trends  for  a  0  to  90  degree  move  of  link  2  (Fig.  9).  For  link  1, 

Pi  and  P2  remain  at  equal  pressures  until  Oi  became  greater  than  35  degrees 

(Fig.  10).   This  is  believed  to  be  caused  by  the  driving  force  on  link  1  of  the 

return  springs.  Once  0i  exceeds  35  degrees,  the  force  added  by  the  springs 


32 


TABLE  I 
VALVE  CONSTANT  CALCULATION  RESULTS 


Input 

current 

(mA) 

Flow 

rate 

(in/s) 

delta 
Ps 

(psi) 

delta 

Pr 

(psi) 

Ks 

Kr 

Joint  1 
Run  1 

2.5 

1.43 
1.15 
.855 

65.9 
54.0 
42.0 

64.7 
74.1 
87.6 

70.4 
62.4 
52.8 

70.4 
53.3 
36.5 

Joint  1 
Run  2 

5.0 

1.86 
1.47 
1.07 

67.4 
56.7 
43.6 

65.0 
75.4 
90.3 

45.3 
39.0 
32.4 

46.1 
33.9 
22.5 

Joint  2 
Run  1 

1.25 

0.410 
0.446 
0.541 

40.7 
47.4 
61.0 

101.3 
98.2 
92.7 

51.4 
51.8 
55.4 

32.6 
36.0 
45.0 

Joint  2 
Run  2 

3.75 

0.95 
1.12 

1.36 

36.0 
40.8 
49.5 

96.3 
90.1 
82.5 

42.2 
46.8 
51.6 

25.8 
31.5 
39.9 
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Figure  9 
Joint  1  Pressure  Plot 
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Figure  10 
Joint  2  PressurePlot 
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decreased  significantly,  and  the  pressure  plots  follow  the  anticipated  pattern. 
The  maximum  change  in  either  Pt  or  P2 ,  for  either  link  was  determined  to  be 
25  psi.  When  incorporated  into  the  valve  model  (Eqn.  15),  the  change  in  the 
flowrate  or  joint  velocity  is  relatively  small,  even  with  a  large  angular, 
change.  These  observations  tend  to  confirm  Lewis'  hypothesis  that  the 
NEPTUNE  is  a  low  performance  robot  arm,  and  as  such,  most  of  the  system 
dynamics  may  be  discounted. 

Because  classical  valve  operation  and  arm  dynamics  were  unable  to  be 
validated  by  actual  component  performance,  it  was  decided  that  an  analytical 
model  of  the  system  would  provide  poor  results.  As  an  alternative,  an 
equivalent  linearization  technique  utilizing  frequency  response  was  selected. 
The  system  was  modeled  empirically  using  magnitude  and  phase  plot. 
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IV.  EMPIRICAL  SYSTEM  CHARACTERIZATION 

A.    EXPERIMENTAL  SETUP 

A  Hewlett-Packard  3562A  Dynamic  Signal  Analyzer  was  used  to  produce 

the  Bode  plots  for  plant  modeling  and  analysis.    This  signal  analyzer  is 

capable  of  scanning  a  desired  frequency  range  with  a  sinusoidal  output,  while 

simultaneously  measuring  and  plotting  transfer  function  magnitude  and 

phase  data.  The  signal  analyzer  output  was  hooked  directly  to  the  servovalve 

electrical  inputs.  The  input  to  the  signal  analyzer  was  taken  from  the  existing 

position  feedback  potentiometers  on  the  manipulator  arm.  A  small  amplitude 

voltage  was  used  to  excite  the  servovalves.     In  this  way  an  equivalent 

linearization  of  the  arm  was  measured. 

Based  on  the  most  likely  operating  regions  of  the  arm,  three  linearization 
configurations  were  selected.  For  link  1,  a  mid-range  point  of  45  degrees  was 
chosen.  It  was  believed  this  position  best  accounted  for  actuator  dynamics, 
link  motion,  and  the  nonlinearity  created  by  the  two  return  springs  on  link  1. 
Three  linearization  points  were  chosen  for  link  2.  They  were  -45,  0,  and  +45 
degrees  (Fig.  11).  These  positions  provided  the  best  coverage  for  link  2 
motion.  Table  II  lists  all  data  acquisition  runs  and  initial  conditions.  Figure 
12  displays  the  anticipated  MIMO  plant  model  to  include  cross  coupled 
dynamics  between  the  joints.  It  was  anticipated  that  each  linearization  point 
would  require  a  different  set  of  transfer  functions. 

Although  the  effects  of  servovalve  fluid  leakage  between  the  spool  and 
housing  were  initially  deemed  negligible  compared  to  the  flow  through  the 
valve,  the  experiments  showed  that  drift  in  the  position  of  the  link  became 
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Run  No. 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 


TABLE  II 
EQUIVALENT  LINEARIZATION  DATA  RUNS 


Ui 

u2 

01 

02 

Function 

(VDC/DC 

(VDC/DC 

(deg) 

(deg) 

Sampled 

offset) 

offset) 

.3/.1 

fixed 

+  45 

-45 

0i/Ui 

.3/.1 

fixed 

+  45 

0 

©l/Ui 

.3/.I 

fixed 

+  45 

+  45 

0!/Ui 

fixed 

.15/.1 

+  45 

-45 

02/U2 

fixed 

.4/.15 

+  45 

0 

02/U2 

fixed 

.2/.1 

+  45 

+  45 

02/U2 

fixed 

.15/.1 

+  45 

-45 

0l/U2 

fixed 

.4/.15 

+  45 

0 

0i/U2 

fixed 

.2/.1 

+  45 

+  45 

0i/U2 

.3/.1 

fixed 

+  45 

-45 

02/Ui 

.3/.1 

fixed 

+  45 

0 

02/U! 

.3/.1 

fixed 

+  45 

+  45 

02/Ui 

01=  45' 
02  =  -45' 

01=45'                               0i=45 

02=    0                                 02  =  45 

yA,     02 

/      / 

// 

/  /  /  /  y 

1 

2     .                                    3 

Figure  11 
Linearization  configurations 
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Figure  12 
Plant  Block  Diagram 


noticeable  when  the  valve  was  hydraulically  charged,  but  electrically 
inactive.  To  negate  this  effect,  the  unactivated  link  fluid  supply  hoses  were 
disconnected  and  plugged  at  the  valve  end.  The  result  of  this  was  that  the 
hydraulic  characteristics  (bulk  modulus  of  elasticity,  fluid  properties, 
pressure  wave  reverberations,  etc.)  of  the  actuating  cylinders  and  supply 
hoses  on  the  unactivated  link  was  included  in  the  equivalent  linearization, 
with  the  servovalve  assumed  to  have  no  leakage. 
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B.     RESULTS 

1.       Link  1 

Figure  13  shows  the  results  of  the  frequency  response  plots  at  0i  = 
+  45  degrees,  02=  -45  degrees.  In  this  run,  02  was  held  constant  and  the 
sinusoid  signal  was  input  to  servovalve  1.  The  output  from  the  joint  1  feedback 
potentiometer  was  compared  to  the  valve  input  signal  in  order  to  obtain  the 
Bode  plots.  The  frequency  range  of  measured  output  was  from  .1-10  Hz.  Curve 
approximations  to  data  (dashed  line)  showed  the  best  fit  in  both  magnitude 
and  phase  plots  was  a  two  pole  system.  One  pole  was  located  at  0  Hz  and  the 
second  pole  was  rooted  at  3.6  Hz.  The  gain  Kp  was  1.1.  The  system  displayed 
an  apparent  transportation  lag  of  0.04  seconds.  This  feature  was  anticipated 
since  all  real  systems,  particularly  mechanical  ones,  have  some  finite  time 
delay  between  signal  input  and  noticeable  output.  The  transportation  lag 
was  constant  for  all  data  runs  of  both  valves,  mainly  because  both  valve  and 
cylinder  assemblies  are  essentially  identical  with  the  supply  hoses  of  about 
the  same  length.  A  trial  first  order  curve  fit  of  the  same  data  with  a  single 
pole  located  at  0  Hz  corresponding  to  a  transfer  function  of  (K/s)  was  also 
tried.  Although  the  approximation  was  relatively  good  at  lower  frequencies, 
higher  frequencies  showed  poor  correlation  in  both  magnitude  and  phase. 
Adjustments  for  transportation  lag  provided  no  significant  improvement  and 
the  overall  result  was  considered  unacceptable.  After  careful  analysis,  the 
following  transfer  function  was  determined  to  most  closely  model  link  1: 


e.      i.i  e-0-04*- 


I 


U  s{s  +  3.6) 


(Eqn  22) 
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Figure  13 
Configuration  1  Linearization  -  0i/Ui  data 
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Linearization  of  link  1  at  02=  0  degrees  (Fig.  14),  and  +45  degrees  (Fig. 
15)  produced  similar  results  as  described  above  with  the  transfer  functions 
being : 


(Eqn.  23) 


el 

lle-O04. 

ui  ~ 

sis  +  3.6) 

el 

1.0  e-0M' 

ux~ 

s(s+  3.6) 

(Eqn.  24) 


respectively.  Since  all  three  transfer  functions  were  within  an  experimental 
errror  factor  of  10%,  it  was  decided  that  one  common  linearization  point  would 
suffice  (Eqn.  22).  Converting  0i  from  volts/sec  to  degrees/sec  gave  the  final 
form  that  would  be  implemented  for  the  link  1  transfer  function. 


©t         100.16  e~0Ms 

—  -  (Eqn.25) 

Ul         s(s  +  327.7) 


2.       Link  2 

Figure  16  shows  data  acquired  from  link  2  with  Oi  fixed  at  +45 
degrees,  and  02  linearized  around  -45  degrees.  Analysis  of  the  plots  revealed  a 
pole  at  0  Hz,  similar  to  link  1  results.  However,  no  signs  of  any  higher 
frequency  poles  were  evident.  A  good  curve  fit  to  data  was  accomplished  with 
a  transfer  function  of  the  form  (K/s).  Once  again  a  transportation  lag  of  0.04 
seconds  was  evident  in  the  phase  plot,  and  the  gain  Kp  was  determined  to  be 
0.21: 
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Configuration  2  Linearization  -  ©i/Ui  data 
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Configuration  3  Linearization  -  Qi/Ui  data 
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Configuration  1  Linearization  -  02AJ2  data 
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©2        .21  e"004' 

—  =  (Eqn.  26) 

u2 


Figure  17  shows  the  results  of  linearization  around  02  =  0  degrees. 
Comparable  results  were  obtained  from  this  data  run,  with  the  transfer 
function  determined  to  be  : 


©2        .26  e-0Ms 

—  =  (Eqn.  27) 

U2  s 


Linearization  at  02  =  +45  degrees  could  not  be  accomplished  using  the 
signal  analyzer  because  of  the  low  DC  signal  offset  required  to  keep  the  link 
in  gravitational  equilibrium.  Limited  data  points  were  taken  manually  with 
Fig.  18  showing  an  approximate  curve  fit.  Although  this  data  was  unusable  in 
providing  an  accurate  model,  it  does  show  correlation  to  the  data  of  the 
previous  two  linearization  points.  The  approximate  transfer  function  for  link 
2  was: 


02  #23e -0.04s 

—  =  ' (Eqn.  28) 

u2 


As  in  the  case  of  link  1,  since  the  results  were  so  closely  grouped,  it  was 
felt  that  only  one  linearization  equation  would  be  needed  to  adequately  model 
the  range  of  the  three  original  set  points.  Conversion  from  volts/sec  to  degree/ 
sec  produced  the  final  form  for  the  link  2  transfer  function  (Eqn.  29). 
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Configuration  2  Linearization  -  02/U2  data 
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©2       {22.5  e"004'} 

—  =  (Eqn.  29) 

u2 

3.  Cross  Joint  Coupling 

In  most  complex  mechanisms,  dynamic  cross  coupling  between  links 

is  a  very  real  factor  that  will  inevitably  increase  the  system  complexity.  In  an 
effort  to  account  for  this  the  plant  model  presented  in  Figure  12  contains  two 
additional  transfer  functions,  0i/U2  and  O2/U1.  These  transfer  functions 
contain  the  cross-coupling  terms  that  generate  an  output  motion  in  link  2 
when  joint  1  is  actuated,  and  vice  versa.  Data  taken  in  runs  7  through  12 
were  used  to  determine  the  cross-coupled  effects.  In  all  runs,  random  data 
scatter  occured  between  -35  and  -70  dB  with  no  correlation  or  consistency.  The 
conclusion  is  that  the  low  velocities  and  accelerations  generated  by  this  arm 
are  not  enough  to  produce  any  significant  cross-coupling  motion.  In  effect  each 
individual  link  can  be  considered  independent.  This  result  is  specific  to  the 
NEPTUNE  arm  performance  and  may  not  be  the  case  for  higher  performance 
systems.  Therefore,  care  must  be  taken  to  test  for  the  existence  of  dynamic 
cross-coupling  before  attempting  similar  simplification. 

4.  Revised  Plant  Model 

With  the  elimination  of  cross-coupling  and  the  decrease  in  the 

number  of  linearization  points,  a  revised  plant  model  is  shown  in  Figure  19. 
The  conversion  of  these  transfer  functions  from  the  Laplace  domain  to  the 
time  domain  yielded : 

9    =  -327.7  Q{  4-  100.2  U{  (t-  0.04)  (Eqn.  30) 

02=  22.5  U2(t-  0.04)  {Eqn.  31) 
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The  acceleration  term  in  Equation  30  strongly  suggests  that  although  the  U2 
input  has  no  cross-coupling  effect  to  0i,  the  inertia  generated  by  the  mass 
translation  of  link  2  is  still  important  to  the  motion  of  link  1.  The  impact  of 
Equation  31  is  that  the  velocity  of  link  2  is  directly  controllable  by  the  input 
voltage  U2,  to  the  servovalve. 


Ui 

100.2    e-0.04s 
.s(s+  327.7) 

Oi 
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22.5e-0.04s 
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02 

Figure  19 
Revised  Plant  Model 
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V.  OPTIMAL  CONTROL  THEORY 

A.    CONTINUOUS  TIME  LINEAR  QUADRATIC  REGULATOR 

The  Linear  Quadratic  Regulator  (LQR)  is  considered  one  of  the  most 
important  and  most  used  optimal  controllers.  The  reason  for  this  is  based  on 
the  fact  that  although  the  derivation  of  the  control  law  is  a  laborious  and 
painstaking  procedure,  the  implementation  is  an  extremely  simple  and 
straightforward  method  to  achieve  multiple  coordinated  inputs.  LQR  control 
assumes  that  the  plant  can  be  modeled  linearly  in  state  space  with: 

x  =  Ax  +  Bu  lEqn.  32) 

Additionally,  the  performance  index  J  must  be  of  the  quadratic  form: 

._/=  /J{\(xrQx)  -  (uTRu)}dt  (Eqn  4) 

As  discussed  earlier,  the  weighting  matrices  Q  and  R  have  special 
definitions  with  respect  to  state  error  and  control  cost.  Further  restrictions 
apply  to  each  weighting  matrix  :  Q  must  be  either  a  positive  semi-definite  or 
positive  definite,  real  symmetric  matrix;  R  on  the  other  hand  is  required  to  be 
a  positive  definite,  real  symmetric  matrix.  Once  Q  and  R  are  reduced  to  the 
proper  form,  a  gain  matrix  K  can  be  found  in  a  two  step  procedure.  The  first 
step  necessiates  the  solving  of  the  reduced  matrix  Riccati  equation  for  the 
unique  positive  definite,  real  symmetric  matrix  P. 

ArP  +  P  A  -  PB  R_1BTP  +  Q  =  0  i.Eqn.33) 

The  matrix  K  will  then  follow  directly  as: 
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K=  R_lB'  I1 


(Eqn.  34) 


With  K  determined,  the  governing  control  law  is  given  as 


u  =  -Kx 


iEqn.5) 


Figure  20  shows  the  corresponding  system  block  diagram. 
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Figure  20 
Optimal  Feedback  Control  Block  Diagram 


In  this  work  the  state  selections  were:  xi  =  Oi  t  X2  =  ©i  ,  X3-O2.  Since 
NEPTUNE  has  only  position  feedback  potentiometers,  joint  1  velocities  (0i) 
were  calculated  using  a  backward  difference  method. 

B.    DISCRETE  TIME  LINEAR  QUADRATIC  REGULATOR 

Because  the  system  controller  was  to  be  implemented  in  software  using  a 
computer,  the  continuous  time  LQR  had  to  be  converted  to  a  discrete  time 
format.   This  required  discretizing  the  plant  model,  the  performance  index, 
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and  the  matrix  Riccati  equation.  This  work  has  been  done  previously  [Refs. 
4,15,16]  and  the  results  are  presented  in  Equations  35  (plant),  38 
(performance  index),  and  39  (gain  matrix)  for  comparison  to  the  continuous 
time  LQR  controller  equations. 

x(n  +■  I)  =  <Px<n)  +  I'u(n)  (Eqn.  35) 


<t>  =  eXh  (Eqn.  36) 


r  ={jh0eAt  dt)  B  {Eqn.  37) 


J=yX      {.\(xr(n)Qx(n))~  (uT  ( n)  R\i)(nj\  dt  (Eqn.  38) 

*—  n  =0 


K  =  (R+  rTPrrl  rrP<D  (Eqn.  39) 


C.     EFFECTS  OF  TRANSPORTATION  LAG 

Transportation  lag  in  the  discrete  time  system  is  an  important  factor  from 

the  viewpoint  of  sample  rate  and  system  stability.   If  the  sample  period  (h)  is 

significantly  shorter  than  the  time  delay  (x),  several  control  periods  may  go  by 

before  there  is  a  corresponding  change  in  the  output.    One  method  to  rectify 

this  problem  is  provided  for  in  Reference  18;  It  requires  the  addition  of  an 

extra  state  for  each  overlapped  sample  period.  With  t  less  than  the  sample 

period,  Equations  40,41,42,43     detail  the     procedure  and  contents  of  the 

modified   *  and  r  matricies.   Actual  determination  of  the  sample  period  that 

was  used  will  be  discussed  in  a  later  section. 
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<t>  =  e 


\h 


{Eqn.  40^1 


r.  =e 


_     A  (A  — i)  /  ri    A* 


{/     "rfOB 


(Eqn.  41) 


r0={/l0*"W  eXtdt\B 


(Eqn.  42) 


x  (n+1) 
u(n) 


<t>  r. 


o  o 


x(n) 
u  (n-1) 


I 


u(n)  (Eqn.  43) 


D.    DISCRETE  PLANT  DETERMINATION 

From  Equations  30,31  the  continuous  time  linear  plant  model  was 

determined  to  be: 


x  = 


1  0  0 

0     -327.7.  0 

0  0  0 


X   + 


0 

0 

100.2 

0 

0    • 

22.5 

U  (Eqn.  44) 


Discretization  of  this  for  a  sample  period  of  0.05  seconds  leads  to: 


X(n-i-l)  = 


1  3.05(10-3)  o 
0  7.66(10-8)  0 
0  0  1 


X(n)  + 


1.44(10-2)    o 
3.03(10-1)   o 
0        1.125 


U  (n)        (Eqn.  45) 


Inclusion  of  the  transportation  delay  produced  the  following  *  and  r  matrices: 


*  = 


1  3.05(10-3)  0      2.03(10-2)  0 

0  7.66(10-8)  0       1.15(10-2)  0 

0            0  10  0.9 

0            0  0           0  0 

0            0  0           0  0 


(  Eqn.  46) 
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r    = 


2.16(10*i  0 

2.94U0-H  0 

0  2.25(10-1) 

1  0 
0  1 


(  Eqn.47) 


E.     DETERMINATION  OF  THE  WEIGHTING  MATRICES 

The  selection  of  the   R  matrix  was  driven  by  the  goal  to  minimize 

hydraulic  power.  To  fit  the  form  of  the  quadratic  performance  index  and  also 

to  force  R  to  be  positive  definite  it  was  decided  to  minimize  Power 2 . 


Power.2  -    P  q.)2 


(Eqn.  48) 


«2     9 
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Joint  2  fit  the  form  exactly  with  ©9=  22.5  ISo-  Because  of  the  second  order 
term  in  the  linearized  equation  for  joint  1,  however,  a  similar  form  did  not 
exist  for  Oi  andUi.  To  alleviate  this,  a  reversion  back  to  the  first  order 
approximation  of  joint  1  motion  was  required  for  the  purpose  of  merit  function 
evaluation  only.  As  discussed  earlier,  a  satisfactory  fit  at  lower  frequencies 
was  found  for  the  (K;s)  joint  1  model,  but  increasing  deviation  at  higher 
frequencies  was  observed.  Actual  observation  of  the  link  in  operation 
indicated  that  the  higher  frequency  dynamics  occured  mainly  in  the 
transition  from  static  equilibrium  to  sustained  motion,  and  vice  versa.  For  the 
majority  of  the  operating  region  much  lower  velocities  were  encountered. 
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Based  on  this  it  was  felt  that,  although  the  approximation  was  not  accurate 
over  the  entire  operating  spectrum,  reasonable  results  could  be  obtained  for 
power  comparisons  with  the  (K/s)  model  for  link  1.  This  led  to  the  following  R 
matrix  based  on  first  order  models  for  both  joints: 


R  = 


80473       0 
0       78400 


(Eqn.  49) 


The  Q  state  weighting  matrix  took  the  form  of  a  reduced  identity  matrix 
I.  The  states  that  were  considered  important  were  0i,  0i,  and  02-  Minimizing 
the  feedback  states  of  past  inputs,  or  Ui(n-l),  was  felt  to  have  no  real  utility 


Q  = 


10     0     0  0 

0      10     0  0 

0      0      10  0 

0      0     0     0  0 

0     0     0     0  0 


(Eqn.  50) 


and  hence  was  ignored.  Although  discussed  only  briefly  in  most  control  texts, 
the  proper  selection  of  Q  is  the  major  factor  in  determining  which  states  the 
controller  will  emphasize  in  optimizing  plant  performance. 

In  order  to  ascertain  the  affect  that  various  magnitudes  of  X  (Eqn.  38) 
have  on  power  minimization,  three  trial  values  were  chosen.  The  first  value  of 
30,  was  picked  specifically  to  show  a  situation  where  both  valve  1  and  2  intial 
voltages  were  in  the  saturated  region.  The  second  value  of  15  placed  the  intial 
voltages  at  the  saturation  boundary.  The  final  value  of  1  produced  valve 
inputs  that  remained  in  the  acceptable  performance  range  throughout  the 
move. 
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F.     CALCULATION  OF  THE  GAIN  MATRIX  K 

With  «,  T,  R,  and  Q  determined,  the  Riccati  equation  was  solved  for  K. 
The  following  values  were  obtained. 


K 


\=i 


3.52U0-3)     1  08(10  ~5 )  °  7.17(10      )  0 

0  0  3.56(10"3)  0  32U0-3) 


iEqn.  51) 


K 


\=15 


1.37U0"2)     4.17(10~5)  °  2.78(10      )  0 

0  0  1.37U0"2)  0  1.24U0-2 


(Eqn.  52) 


K 


\  =  30 


1  93(10~2)    5.89U0-5)  °  3.94(10-4)  0 

0  0  1.93(10"2)  0  1.74(10"2) 


(Eqn.  53) 
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VI.    MICRO-COMPUTER  IMPLEMENTATION 

A.    CONTROLLER  SAMPLE  RATE 

Determination  of  the  correct  controller  sample  rate  for  the  IBM  PC  AT  in 
conjunction  with  the  NEPTUNE  arm  was  a  two  part  process.  First,  an  offline 
software  test  was  conducted  to  determine  the  time  required  for  the  execution 
of  the  complete  control  algorithm.  After  this  was  accomplished,  an 
examination  of  the  arm  model  was  performed  to  predict  the  minimum  sample 
rate  that  was  needed  perform  the  ADC,  control,  and  DAC  while  still 
maintaining  system  stability. 

The  offline  software  test  consisted  of  coding  the  entire  control  loop  in 
Fortran,  and  compiling.  This  test  code  was  then  iterated  5000  times,  and  the 
elapsed  time  was  recorded.  The  results  of  this  test  produced  an  average 
computer  control  cycle  period  of  2.57  milli-seconds.  As  a  general  rule  of 
thumb,  it  was  decided  that  the  control  computation  time  should  not  exceed 
60%  of  the  computer  control  sample  period.  The  reasoning  behind  this  was 
that  system  overhead  time  would  be  required  for  sensor  acquisition, 
conversion,  and  computer  "housekeeping"  chores.  For  LABPAC  this  overhead 
time  is  approximately  0.8  milli-seconds,  making  the  minimum  computer 
sample  period  3.37  milli-seconds. 

Shannon  sampling  criteria  was  used  to  determine  the  minimum  sample 
frequency  required  to  enable  reconstruction  of  the  sampled  motion  [Ref.  1). 
Since  the  highest  frequency  pole  was  associated  with  joint  1  at  3.6  Hz,  the 
lowest  acceptable  sample  frequency  for  the  system  was  set  at  7.2  Hz,  with  a 
recommended  sample  frequency  of  36  Hz.  This  equated  to  a  sample  period  for 
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system  stability  of  between  .028  (recommended)  and  .139  (required)  seconds. 
Use  of  the  recommended  sample  rate,  however,  meant  that  the  transportation 
lag  of  0.04  seconds  would  extend  over  two  sample  periods.  A  compromise  .05 
second  (50  msec)  sample  period  (20  Hz)  was  implemented,  thus  making  the 
duty  cycle  (3.37  msec)  less  than  10%  of  the  sample  rate. 

B.    POWER  DATA  COMPARISONS 

A  total  of  four  comparison  runs  were  performed  with  the  arm.  The  first 
three  runs  used  the  gain  matricies  that  corresponded  to  the  three  previously 
discussed  values  of  A.  The  last  run  utilized  the  solenoid  control  as  originally 
configured  on  the  NEPTUNE  to  provide  a  baseline  for  comparison.  The  start 
point  for  each  run  was  0i  =  2O  degrees,  02  =  -60  degrees;  with  the  end  point 
being  0i  =  70  degrees,  02  =  30  degrees.  Figures  21,22  show  time-motion  plots 
of  joints  1  and  2.  Each  link  shows  expected  results,  with  the  optimal  control 
runs  approaching  their  desired  set  points  assymtotically.  It  should  be  noted 
here  that  in  the  case  of  A  =1,  the  target  position  was  never  actually  reached. 
02  finally  ended  3  degrees  short  of  the  desired  position,  because  vaive  input  - 
Kx  was  in  the  deadband  voltage  for  that  small  error  condition  .  In 
comparison,  the  solenoid  runs  approached  the  target  positions  with  constant 
velocity,  and  motion  was  slowed  by  use  of  the  fluid  restrictors.  When  the  link 
was  within  the  error  deadband,  motion  was  stopped  by  shutting  the  solenoid 
valve. 

Figures  23,24,25  compare  total  hydraulic  power  consumption  for  the  arm 
(  E(Pii-P2i)qi  )  with  solenoid  power  consuption,  for  the  specified  \'s.  The 
negative  power  values  that  appear  for  the  solenoid  and  \=1  move,  are 
attributable  to  the  fact  that  the  joint  1  return  springs  were  providing  all  forces 
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Figure  21 
Joint  1  position  plot 
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Figure  22 
Joint  2  position  plot 
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Figure  23 
Total  arm  power  -  X  =  1 
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Figure  24 
Total  arm  power  -  A.=  15 
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Figure  25 
Total  arm  power  -  A.  =  30 
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necessary  for  motion  of  link  1  until  0i  exceeded  35  degrees.  This  phenomenon 
was  not  evident  for  A  =  15,  or  A  =  30.  As  A  increased  from  1  to  30,  the 
maximum  instantaneous  power  used  by  the  arm  was  correspondingly  higher 
reflecting  increased  weight  on  error.  Accurate  values  for  the  actual  power 
used  by  the  arm  were  difficult  to  determine,  largely  because  of  the  noise  from 
the  position  feedback  potentiometers,  and  random  errors  in  ADC.  However,  a 
general  trend  can  be  discerned  that  shows  A  =  1  had  a  lower  maximum  power 
requirement  than  did  the  solenoid.  The  sawtooth  like  curve  of  the  solenoid 
operation  correlates  with  fluid  restrictor  closing.  The  power  requirements  for 
all  optimal  control  moves  appear  to  transition  smoother,  from  high  to  low 
values.  This  can  be  advantageous,  especially  in  hydraulic  systems,  where 
sudden  shocks  are  not  desired. 

Total  system  power  plots  (power  taken  from  accumulator  with  Ps 
assumed  constant)  show  results  consistent  with  arm  power  results  .  A  =  30  and 
15  plots  (Fig.  26,27)  show  that  the  power  used  by  the  system  always  exceeds 
solenoid  power  requirements.  In  contrast,  the  power  requirements  for  A  =  1 
somewhat  closely  matches  that  of  the  solenoid  (Fig.  28).  The  major  deviation 
is  the  initial  power  surge  for  the  servovalve  operation.  Once  again, 
transitions  in  power  levels  were  smoother  in  optimal  control  than  solenoid 
control. 
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Figure  26 
Total  system  power  -  X  =  30 
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Figure  27 
Total  system  power  -  A.  =  15 
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Figure  28 
Total  system  power  -  \  =  1 
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vn.  CONCLUSIONS 

A.  SYSTEM  LINEARITY 

The  intial  evaluation  of  the  equations  of  motion  for  the  NEPTUNE 
system  were  based  on  the  assumption  that  the  rotational  inertia  of  the  arm 
structure  was  dominant.  This  in  fact  was  not  the  case  with  the  NEPTUNE. 
With  the  servovalves  installed,  the  dominant  mechanism  for  determining 
motion,  was  the  combination  servovalve  and  actuator  assembly.  Even  at 
relatively  low  pressures,  the  NEPTUNE  behaves  linearly.  In  essence,  the 
hydraulic  system  has  enough  power  to  override  much  of  the  low  performance 
dynamics  that  might  have  been  expected. 

B.  SERVOVALVE  PERFORMANCE 

The  data  that  was  obtained  for  the  servovalve  did  not  reflect  the 
anticipated  classical  valve  model.  The  reason  for  this  is  believed  to  be  the 
relatively  low  system  pressure  and  the  distance  between  the  servovalves  and 
pressure  sensors.  It  was  felt  that  by  correcting  these  factors,  better 
correlation  may  be  obtained.  Nevertheless,  equivalent  linearization 
techniques  were  used  to  adequately  model  both  the  valve  and  arm 
performance  as  configured  in  the  present  study. 

C.  MICRO-COMPUTER  OPTIMAL  CONTROL 

Micro-computer  optimal  control  is  a  simple  to  implement,  feasible  control 
scheme.  The  control  system  can  be  "tuned"  to  the  desired  performance  level 
by  the  proper  choice  of  weighting  matrices  Q  and  R.  Additionally,  a  smooth 
transition  in  motion  and  hydraulic  power  requirements  can  be  expected. 
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D.    SAMPLE  RATE 

Since  the  NEPTUNE  was  of  such  low  performance,  the  maximum  sample 
rate  available  with  the  IBM  PC-AT  was  adequate.  The  driving  concern  in 
sample  rate  determination  was  the  proper  handling  of  the  transportation 
delay.  Finally,  although  the  selected  sample  rate  was  less  than  a 
recommended  five  times  the  Shannon  sampling  criteria,  system  stability  was 
not  affected. 
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Vm.  RECOMMENDATIONS 

Continue  the  use  of  NEPTUNE  plus  IBM  PC-AT  as  a  test  bed  for 
controller  development.  Three  future  projects  are  recommended  as  follow  on 
work.  First,  a  three  dimensional  optimal  controller  should  be  implemented  in 
order  to  utilize  the  full  range  of  motion  on  the  NEPTUNE  arm.  Second, 
comparison  studies  of  optimal  versus  Proportional-Integral-Derivative  (PID) 
control  should  be  conducted,  to  provide  a  better  baseline  evaluation  against  a 
common  type  of  controller  found  in  industry.  Finally,  the  NEPTUNE  and 
IBM  PC-AT  should  be  configured  to  perform  more  intelligent  types  of  control 
such  as:  adaptive,  obstacle  avoidance,  etc. 

Additionally,  the  following  hardware  improvements  should  be  made  to 
the  system.  First,  it  is  recommended  that  the  system  pressure  be  increased  to 
that  specified  by  the  valve  manufacturer,  and  the  system  should  be 
remodeled.  It  is  also  recommended  that  the  rack  and  pinion  joint 
transmission  gears  be  redesigned  to  remove  or  minimize  backlash.  Lastly, 
digital  shaft  encoders  should  be  installed  to  increase  the  accuracy  of  joint 
position  feedback. 
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APPENDIX 
OPTIMAL  CONTROL  PROGRAM 

$include:'labhead.for' 

Sincludei'master.for' 

$include:'dadio.for' 

dimension  chan(6),databufT6,2500) 

real*4  volt,runt,terml,term2,svltlr,svlt2r,kopt(2,5),deltjl,deltj2 

real*4stptl,stpt2,endptl,endpt2,  lambda,  thdotl 

logical  flagl,flag2 

c  constants 

thlzero  =  41 
th2zero  =  1300 
thlnint  =  2620 
th2nint  =  3835 
thlrng  =  thlnint-thlzero 
th2rng  =  th2nint-th2zero 
count  =1 
runt  =  40 
nsweep  =  runt*20 
dither  =  10 

c  set  channel  numbers  to  be  scanned 

chan(l)  =  l 
chan(2)  =  2 
chan(3)  =  6 
chan(4)  =  7 
chan(5)  =  8 
chan(6)  =  9 


c  reset  and  initialize  labpac 


result  =  labpac(RESET) 

result  =  labpac(ATOD,12,AlINIT) 

c  initialize  dadio  board  for  digital  to  analog  output 

result  =  labpac(DA,4,AOINIT) 
result  =  labpac(7,125,TIMER,SWINIT) 
result  =  labpac(0,0,AORAW) 
result  =  labpac(  1 ,0,AORAW) 


c  read  in  the  starting  and  ending  joint  positions 
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400       write  (*,*)  'ENTER  JOINT  1  STARTING  POSITION  (DEGREES) :' 
read  (V)stptl 

write  (V)  "ENTER  JOINT  1  ENDING  POSITION  (DEGREES) :' 
read  (*,*)  endptl 

write  (*,*)  'ENTER  JOINT  2  STARTING  POSITION  (DEGREES) :' 
read  (*  *)  stpt2 

write  (*,*)  'ENTER  JOINT  2  ENDING  POSITION  (DEGREES) :' 
read  (*,*)  endpt2 
write  (*,*) 
write  (*,*) 

c  read  in  the  array  of  control  feedback  constants 

write(V)  'ENTER  THE  FILE  NAME  THAT  CONTAINS  THE  K 
+  MATRIX' 
write  (*,*) 

open  (12,  file  = '  '.status  = 'old') 
read(12,*)  lambda 
do  50  i  =  1,2 
50        read  (12,*)  kopt(i,l),kopt(i,2),kopt(i,3),kopt(i,4),kopt(i,5) 


c  move  joint  1  and  2  to  starting  position 

pause  'TURN  ON  MASTER  SERVO  POWER  SWITCHES' 
flagl  =  .  false. 
flag2  =  . false.  ■ 
term2  =  stptl/90.*thlrng 
iterm2  =  ifix(term2) 
100       thl  =  labpac(l,AIRAW) 
delthl  =  thl-iterm2 
svitlr  =  -(delthl*2047.)/thlrng 
svoltl  =  ifix(svltlr)  +  dither 
if(svoltl  .It. -2047)  then 

svoltl  =-2047 
endif 
if  (svoltl  .gt.  2047)  then 

svoltl  =  2047 
endif 
if  (abs(delthl)  .le.  15)  then 

result  =  labpac(0,0,AORAW) 

flagl  =  .true. 
else 

result  =  labpac(0,svoltl,AORAW) 

flagl  =  .false, 
endif 

th2  =  labpac(2,AIRAW) 
terml  =  stpt2*(th2rng/90.) 
iterml  =  inx(  terml) 
delth2  =  th2-th2nint-iterml  +  thl 
svlt2r  =  -(delth2*2047.)/th2rng 
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svolt2  =  ifix(svlt2r)  +  dither 
if(svolt2  .It.  -2047)  then 

svolt2  =  -2047 
endif 
if  (svolt2  .gt.  2047)  then 

svolt2  =  2047 
endif 
if  (abs(delth2)  .le.  15)  then 

result  =  labpac(  1 ,0,AORAW) 

flag2  =  .true. 
else 

result  =  labpac(  1  ,svolt2,  AOR  AW) 

flag2  =  .false, 
endif 
if  (flagl  .and.  flag2)  then 

goto  200 
endif 

dither  = -dither 
goto  100 


c  commence  optimal  control  and  data  acquisition  run 

200      write(V) 

pause  'READY  TO  COMMENCE  OPTIMAL  CONTROL  MOVE' 

c  convert  theta  2  from  global  to  relative 
endpt2  =  180.-endptl  +endpt2 

swleft  =  nsweep- 1 

lstoutl=0. 

Istout2  =  0. 

result  =  labpac(l,NSWEEP,6,CHAN,DATABUF,AISWST) 


c  optimal  control  algorithm  with  20th  of  a  second  time  delay 
c  between  control  sweeps 

300       result  =  labpac(SWLEFT,AISTAT) 

deltjl  =  float(databufU,count)-thlzero)*90./float(thlrng)-endptl 
deltj2  =  float(databufT2,count)-th2zero)*90./float(th2rng)  +  90. 
+  -endpt2 

ift((abs(deltjl).le.  1.)  .and.  (abs(deltj2).le.  1.))  .or.  (swleft 
+  .eq.  0))then 

result  =  labpac(0,0,AORAW) 
result  =  labpac(  1 ,0,  AORAW) 
result  =  labpac(AISWAB) 
stop 
endif 

thdotl=float(databufTl,count-l)-databuf\l,count))*1800./float 
+  (thlrng) 
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svltlr  =  -kopt(l,l)*deltjl-kopt(l,2)*thdotl-kopt(l,3)*deltj2- 
+  kopt(l,4)*lstoutl-kopt(l,5)*lstout2 

svoltl  =  ifix(svltlr*2047/10.  +  dither) 

svlt2r  =  -kopt(2,l)*deltjl-kopt(2,2)*thdotl-kopt(2,3)*deltj2- 
+  kopt(2,4)*lstoutl-kopt(2,5)*lstout2 

svolt2  =  ifix(svlt2r*2047/10.  + dither) 

result  =  labpac(0,svoltl,AO  RAW) 

result  =  labpac(l,svolt2,AORAW) 

lstoutl=svltlr 

lstout2  =  svlt2r 

count  =  count +  1 

swleft  =  swleft-l 

dither  = -dither 

go  to  300 

end 
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LIST  OF  VARIABLES 

A     The  cross-sectional  area  (in2) 

c     The  distance  from  the  neutral  axis  to  the  most 
extreme  fiber  (in) 

~d~     Support  tube  mean  diameter  (0.9  in) 

E     Young's  Modulus  of  Elasticity  (29,000,000   lbs/in2) 

Fr  ,  Fr  ,  Fr    Resultant  forces  for  an  arbitrary  point 
x    Y  onzthe  support  tube  cylinder  (lbs) 

F1'F2'F3       Load  forces  at  the  penetrator  tip  (lbs) 

G     Shear  Modulus  of  Elasticity  (lbs/in) 

I     The  moment  of  inertia  of  the  cross-section  with 
respect  to  the  neutral  axis  (in4) 

J     The  polar  moment  of  inertia  of  the  cross-section 
relative  to  its  centroid  (in4) 

L21'L22  The  distance  from  the  centerline  of  the  drill  to 
the  critical  points  on  the  support  cylinder  surface 
(L2i  =  2.75  in  and  L22  =  6.75  in) 

M     Bending  moment  (in-lbs) 

Mr  ,Mr  ,  Mr     The  resultant  moments  for  an  arbitrary  point 
x    y     onzthe  cylinder  (in-lbs) 

Mi,M2,M3    The  load  moments  at  the  penetrator  tip  (in-lbs) 

P     Transmitted  force  (lbs) 

Q     The  first  moment  of  the  area  with  respect  to  the 
neutral  axis  (in^) 

r  Support  tube  mean  radius  (0.45  in) 

r^  Inside  radius  of  the  support  cylinder  (0.4  in) 

r0  Outside  radius  of  the  support  cylinder  (0.5  in) 

t  Radial  thickness  of  the  support  cylinder  (0.1  in) 
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T     Torque  (in-lbs) 

V     Transverse  shear  force  (lbs) 

&     The  angle  between  a  force  vector  and  the  x  axis 
(degrees) 

3     The  angle  between  a  force  vector  and  the  y  axis 
(degrees) 

f     Shear  strain  (in/in) 

*xv'^YZ'  zx  T*ie  shear  strain  where  the  first  subscript 
denotes  the  plane  face  on  which  the  strain  acts  and 
the  second  the  direction  on  that  face  (in/ in) 

^yz  (L2l)      T^e  shear  strain  at  the  point  B  a  distance  L«2x 
B     from  the  centerline  of  the  penetrator  (in/ in) 

6     The  angle  between  a  force  vector  and  the  Z  axis 
(degrees) 

t,  Normal  strain  (in/in) 

fiX'£y£z  T^e  normal  strain  where  the  subscript  denotes 
the  plane  face  on  which  the  strain  acts  (in/in) 

£yZ  A  strain  measured  by  a  rosette  gage  element  in  the 
Y-Z  plane  (in/m) 

€-v^(L2]_)  The  normal  strain  at  a  point  A  on  the  cylinder  a 
distance  L2i  from  the  centerline  of  the  penetrator 
(in/ in) 

*y  (L22)     T^e  normal  strain  at  a  point  A  on  the  cylinder  a 
A      distance  L22  from  the  centerline  of  the  penetrator 
(in/in) 

€y  (L21)     The  normal  strain  at  a  point  B  on  the  cylinder  a 
A      distance  L2i  from  the  centerline  of  the  penetrator 
(in/in) 

€y  (L22)     T^e  normal  strain  at  a  point  B  on  the  cylinder  a 
A      distance  L22  from  the  centerline  of  the  penetrator 
(in/ in) 

£y  (L2i)     The  normal  strain  at  a  point  C  on  the  cylinder  a 
B      distance  L21  from  the  centerline  of  the  penetrator 
(in/in) 
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€yz  (L2l)      Tlie  Y~z  Plane  strain  at  a  point  B  on  the 
B     cylinder  a  distance  L21  from  the  centerline  of  the 
penetrator  (in/ in) 

£„  ... 

B(L2l)     The  normal  strain  at  a  point  B  on  the  cylinder  a 
distance  L21  from  the  centerline  of  the  penetrator 
(in/ in) 

0      An  angle  in  the  X-Z  plane  measured  from  the  X  axis 
(degrees) 

V      Poison's  ratio  (0.29) 

<r      Normal  strain  (lbs/in) 

o^,<7y,<J*z      The  normal  stress  where  the  subscript  denotes 
the  plane  face  on  which  the  strain  acts  (lbs/ in) 

^y  ,<TV    ,<rvr   The  norTnal  stress  acting  on  the  Y  plane  face  at 
A   B  the  points  A,  B,and  C 

^(r,^   The  normal  stress  acting  on  the  Y  plane  face  as  a 
function  of  the  variables   r  and  6  (in/in) 

*&      Shear  stress  (lbs/in) 

tyx  ,%rZ    r^yxr        The  shear  stress  at  the  points  A,  B,  and  C 

A    3where  the  first  subscript  denotes  the  plane  face 

on  which   the   stress   acts   and   the  second   the 
direction  on  that  face  (in/ in) 
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I.   INTRODUCTION 

A.  PROBLEM  DESCRIPTION 

An  aircraft  skin  penetrator/applicator  was  developed  by 
the  USAF  as  a  device  for  penetrating  an  aircraft  fuselage 
and  injecting  Halon  into  the  interior  of  the  aircraft  to 
extinguish  fires  [Ref .  1]  .  The  guidelines  for  the  present 
study  and  other  related  Naval  Postgraduate  School  work  were 
to  investigate  the  requirements  for  a  robot  arm  to 
manipulate  such  a  tool.  The  purpose  of  this  work  was  to 
identify  the  loads  that  the  support  structure  would  have  to 
sustain  during  a  drilling  operation  with  the  device.  In 
general,  the  loads  are  functions  of  the  material  character- 
istics' of  the  drilling  surface  and  the  drill  parameters, 
i.e.,  type  drill  bit,  length  of  the  drill,  feed  rate,  and 
angle  of  attack.  For  this  investigation,  the  drill  bit 
parameters  were  predetermined  from  previous  work  at  the 
Naval  Postgraduate  School  [Ref.  2].  Sheet  aluminum  was  used 
to  simulate  the  aircraft  fuselage. 

B.  APPROACH 

The  approach  was  to  experimentally  determine  the  drill 
tip  loads  by  selective  support  structure  strain  measure- 
ments. The  analytical  method  used  was  to  resolve  the  loads 
at  the  drill  tip  into  six  force  and  moment  components  about 
an   orthonormal   coordinate   axes   system  (Figure  1) .   These 

82 


M2(t) 


t 


F3(t) 


F  (t)  A  Mo( 
2 


X 


Figure  1.   Penetrator  Forces  and  Moments 
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loads  create  a  state  of  strain  on  the  outer  surface  of  the 
support  structure  which  was  analyzed  by  conducting  a  static 
force/moment  equilibrium  analysis  and  a  solid  mechanics 
analysis.  The  resulting  strain  equations  were  evaluated  to 
determine  a  set  of  selective  strain  measurements  which  would 
facilitate  segregation  of  the  individual  effects  of  each 
drill  tip  load.  The  analytical  expressions  for  the  strain 
at  those  locations  were  then  used  to  find  analytical 
expressions  for  the  force/moment  loads.  In  this  way,  the 
measured  strains  on  the  support  structure  were  used  to 
estimate  the  six  forces  and  moments  at  the  drill  tip  which 
are  experienced  during  a  drilling  operation. 

C.   PREVIOUS  WORK 

Reference  2  is  a  previous  study  of  the  same  problem  with 
the  assumption  that  the  effects  of  certain  forces  and 
moments  were  negligible  (Figure  2).  Notice  that  only  two 
forces  (F^  and  F3)  and  one  moment  (M3)  were  previously 
assumed  to  be  significant.  The  same  methodology  for 
determining  analytical  expressions  for  the  force/moment 
components  as  functions  of  measured  strain  quantities  was 
used  in  both  studies,  but  the  present  study  removes  this 
restrictive  assumption  on  the  load  magnitudes.  The  previous 
work  also  used  an  Operational  Test  and  Evaluation  drill  to 
make  measurements,  while  this  investigation  used  a  new  drill 
supplied  by  the  NSWC/White  Oak  laboratory. 
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Figure  2.   Penetrator  Forces  and  Moments  of  Previous  Work 
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II.   EXPERIMENTAL  EQUIPMENT 

A.   PENETRATOR  ASSEMBLY  AND  FOUNDATION 

1.  Aircraft  Skin  Penetrator/Applicator 

An  AMETEK,  Offshore  Research  and  Engineering 
Division  model  ASP/AA-1  skin  penetrator/ applicator  was 
utilized  (Figure  3)  .  For  this  work  the  carrier  and  bottle 
assemblies  and  the  Halon  ball  valve  were  removed  to 
facilitate  installation  on  a  milling  machine  bed.  The  drill 
was  designed  to  be  powered  by  compressed  air  carried  by  the 
operator.  In  the  present  work,  a  source  of  compressed  shop 
air  (165  psig)  was  used  to  power  the  drill. 

2 .  Instrumented  Support  Fixture 

A  support  fixture  was  manufactured  from  mild  steel 
which  consisted  of  a  central  cylinder  with  top-  and  bottom 
support  plates  welded  to  its  external  radius.  The 
penetrator  was  attached  to  the  top  plate  with  eight 
self-locking  helicoil  fasteners.  The  bottom  plate  was 
bolted  to  a  milling  machine  bed.  Five  strain  gages  were 
mounted  on  the  cylinder  surface  to  facilitate  the 
measurement  of  the  cylinder  deformation  during  a  drilling 
operation. 

3 .  Milling  Machine 

A  Milwaukee  model  H  milling  machine,  with  mill 
removed,  provided  a  movable  foundation  for  the  penetrator. 
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Two  Superior  Electronic  synchronous/stepping  motors,  72  RPMf 
120  volts  and  60  hertz,  were  installed  on  the  milling 
machine's  transversing  gear  assemblies  and  were  controlled 
by  3 -position  switches.  These  motors  provided  for  movement 
along  the  Y  and  Z  directions  at  a  constant  speed  of  14 
inches/minute . 

B.  SIMULATED  AIRCRAFT  FUSELAGE 

1.  Test  Panels 

As  in  Reference  2,  the  aircraft  fuselage  was 
simulated  by  test  panels  constructed  of  7075-T6  aluminum 
0.050  inch  sheets.  The  panels  were  rolled  to  the  shape  of 
the  support  frame . 

2 .  Support  Frame 

The  support  frame  with  a  test  panel  installed  is 
shown  in  Figure  4 .  The  frame  was  bolted  to  the  floor  and 
leveled.  Due  to  the  manner  in  which  the  support  frame  was 
constructed,  the  minimum  angle  of  attack  achievable,  without 
penetrating  the  support  frame,  was  8.48  degrees  between  the 
drill  axis  and  the  sheet  aluminum. 

C.  STRAIN  MEASUREMENT/RECORDING  EQUIPMENT 
1.   Strain  Gages 

Linear  gages  were  installed  at  critical  points  on 
the  support  structure.  These  gages  were  Micro-Measurement 
Division  series  CEA-06-125UN-350  constantan  strain  gages 
with  a  gage  factor  of  2.065  +/-  0.5%  ,  and  a  resistance  of 
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3  50  ohms  +/-  0.3%.  A  rosettee  strain  gage  was  used  to 
determine  shear  force;  it  was  a  Micro-Measurement  Division  3 
element  4  5  degree  rectangular  stacked  rosette,  series 
WK-06-120WR-350,  superimposed  K-alloy,  strain  gage  with  a 
gage  factor  of  2.065  and  a  resistance  of  350.0  ohms  +/- 
0.4%. 

2 .   Recording  Equipment 

Recording  of  the  strain  gage  outputs  was 
accomplished  with  an  Astro-Med,  Inc.  8  channel  strip  chart 
recorder  via  an  Ectron  amplifier  assembly. 
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III.   ANALYTICAL  METHOD 

A.   BACKGROUND 

1.   System  Model 

The  system  model  was  a  rigid  body  assembly  which 
consisted  of  a  pneumatic  drill  mounted  on  a  rigid,  movable 
support  structure.  The  support  structure  was  a  thin-walled 
tube  (diameter  >>  thickness)  welded  to  support  plates  top 
and  bottom  and  fixed  to  a  rigid  milling  machine  bed. 

2. •  Theoretical  Approach 

The  first  step  was  to  conduct  a  static  analysis  of 
the  system  to  determine  the  force/moment  resultants  at  any 
arbitrary  point  on  the  support  tube  as  a  function  of  the 
drill  tip  loading.  Following  this,  a  solid  mechanic 
analysis  was  conducted  to  determine  the  particular  effect  of 
each  force/moment  component.  Using  the  theory  of 
superposition,  the  six  effects  were  added  to  determine  the 
net  effect  of  the  force/moment  components.  These  were  then 
converted  to  functions  of  drill  tip  loading  by  utilizing 
the  equilibrium  equations  of  the  static  analysis.  The  next 
step  was  to  apply  the  fundamental  triaxial  stress-strain 
relationships.  This  provided  analytical  expressions  for  the 
strain  at  any  arbitrary  point  on  the  tube  as  a  function  of 
the  loads  at  the  drill  tip.  The  final  step  was  to  evaluate 
the   reverse   case;   that   is,   to   determine   analytical 
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expressions  of  the  force/moment  loads  of  the  drill  tip  as 
functions  of  the  strains  at  an  arbitrary  point  on  the  tube. 
3 .   Theory  of  Superposition 

The  theory  of  superposition  for  linear  systems 
states  that  the  net  stresses  of  a  system  loaded  with  a 
combination  of  forces  and  moments  can  be  analyzed  by 
determining  the  resulting  stresses  of  each  load  and  adding 
the  like  stress  components.  Therefore,  an  analysis  of  the 
normal  strains  and  shear  strains  resulting  from  the  reaction 
loads  of  Figure  1  at  an  arbitrary  cut  on  the  tube  was 
analytically  evaluated  for  each  load  and  then  combined  to 
determine  the  net  effect. 

B.   STATIC  SOLID  MECHANICS  ANALYSIS 
1 .   Static  Analysis 

A  free  body  diagram  of  the  system  with  static 
loading  is  presented  in  Figure  5.  The  static  equations  are: 

£FX  =  0:  PX  •+  Fr  =0  ■*  Fr  =  -FX 
J  Fy  =  0:  F2  +  Fr^  =  0  +  Fr^  =  -F2 
IFZ.  -■<>■:  '  Fz  +  Fr  -  0   •»  Frz  =  -F3 

A     A         1     1     A  * 

£Mr  =  0:   (L2j-L1K)X(F1i+F2j+F3k)  +  (Mi+M^)  i 

+  (M2+Mr  )j  +  (M3+Mr  )k  =  0 

y        z 

Equating  the  like  components  equal  to  zero  yields: 
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i  components:   L2F3  +  L1F2  +  Ml  +  Mr*  =  ° 

-  Mrx  =  -(L2F3+L1F2+M1) 

j  components:   ~LiFi  +  M2  +  Mr   =0   ■*  Mr   =  L^F^  -  M2 

k  components:   -L2Fi  +  M3  +  Mr   =0   -*■  Mr   =  L^F]^  -  M3 

z  z 

2 .   Load-Stress  Fundamental  Relationships  TRef.  31 

For  axial  tension,  the  normal  stress  is  given  by 


-  -  I 


where : 

P  =  the  transmitted  force 

A  =  the  cross-sectional  area  (A  =  irtd) 

For  axial  bending,  the  normal  stress  is  given  by: 


•  -!§ 


where: 

M  =  the  bending  moment 

c  =  the  distance  from  the  neutral  axis  to  the 
most  extreme  fiber  (c  =  r0) 

I  =  the  moment  of  inertia  of  the  cross-section 
with  respect  to  the  neutral  axis 

r-r      *t        4     4M 

[I  =  4(ro"ri  )] 

For  transverse  loading,  the  shearing  stress  is  given 
by: 
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r   =  5£ 

1    Jt 


where 


V  =  the  transverse  shear  force 

Q  =  the  first  moment  of  the  area  of  the  semi- 
circular section  with  respect  to  the 
neutral  axis 


[Q=  K3-ri3)l 


t  =  the  radial  thickness  of  the  cross-section 

J  =  the  polar  moment  of  inertia  of  the  cross- 
section  relative  to  its  centroid 


(j  =  Ttd3) 


For   torsion   of   thin-walled   tubes,   the   shearing 
stress  is  given  by: 


*  =  S 


where 


T  =  the  applied  torque 

r  =  mean  radius 

J  =  the  polar  moment  of  inertia  of  the  cross- 
section  relative  to  its  centroid 


(J  =  Jtd3) 
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3 .   Stress  Analysis 

The  state  of  stress  at  a  point  on  the  tube  is 

described  by  nine  values  of  stress.   The  following  is  an 

analysis   of   each   reaction   load   to   determine   their 

contribution  to  these  stress  values. 

For  Mr  ,  Figure  6  depicts  the  loading  condition, 
z 

with  the  resulting  stress  distribution.  Notice  that  the 
stress  at  A  and  C  is  zero,  and  at  B  is  maximum.  These 
critical  points  will  be  important  in  all  the  stress  and 
strain  relationships  which  follow.  The  stress  distribution 
is  given  by: 

M   r  cos  0 

r   o 

<r  (r,6)  =      T 

y  I 

For  Mr  ,  Figure  7  depicts  the  loading  condition, 
with  the  resulting  stress  distribution.  The  stress 
distribution  is  given  by: 

.   M   r  sin  $ 

r   o 

<r(r,e)   =  -2— j 

For  Mr  ,  Figure  8  depicts  the  loading  condition, 
with  the  resulting  stress  distribution.  The  stress 
distribution  is  given  by: 

M   r 

r 

«  =     y 
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Figure  6. 


Mr   Load  and  Stress  Distribution 
z 
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Mr   Load  and  Stress  Distribution 
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Figure  8.   Mr   Load  and  Stress  Distribution 
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For  Fr  ,  Figure  9  depicts  the  loading  condition, 
z 

with   the   resulting   stress   distribution.     The   stress 
distribution  is  given  by: 

F   Q 
rz 

T   -    — - — 

1      Jt 

For  Fr  ,  Figure  10  depicts  the  loading  condition, 
with  the  resulting  stress  distribution.  The  stress 
distribution  is  given  by: 

F   Q 

Jt 

For  Fr  ,  Figure  11  depicts  the  loading  condition, 
with  the  resulting  stress  distribution.  The  stress 
distribution  is  given  by: 

F 
r 

Now,  utilizing  the  theory  of  superposition  the  net 
stress  at  points  A,  B,  and  C  can  be  determined. 
Additionally,  utilizing  the  static  equilibrium  equations  of 
Section  II.B.l,   the  following  six  stress  equations  are 


produced: 

Point  A: 

^A 

M       r         F 
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4 .  Stress-strain   Fundamental   Relationships 

ex     =    I[ax-v(ay  +  az)] 
*y     =     f[ay"  v(Gx  +  az)] 

*2     =     f[az"  V(°x  +  ay)]  ■  - 

*xy     =    VG 

2f       =    r  /g 

yz  y2^ 

•zx  zx 

5.  Determination   of   Strains 

As  stated  in  the  theoretical  approach,  Section 
II. A. 2,  the  next  level  of  analysis  was  to  determine 
analytical  expressions  for  the  strains  at  points  A,  B,  and  C 
as  functions  of  the  loads.    This  was  accomplished  by 
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substituting  the  expressions  for  stresses  of  Section  II. B. 3 
into  the  fundamental  triaxial  stress-strain  relationships  of 
Section  II. B. 4.   Table  I  summarizes  the  results.   Appendix  A 
contains  the  Fortran  computer  code  for  these  equations. 
6.   Determination  of  the  Drill  Tip  Loads 

Figure  12  shows  the  locations  of  the  strain  sites  of 
interest.  Point  A  is  a  point  on  the  outer  surface  of  the 
tube  in  the  Y-Z  plane.  L2j_  is  the  distance  from  the 
centerline  of  the  drill  to  any  arbitrary  Z-X  plane  cut  of 
the  tube  (e.g.,  horizontal  cut)  through  point  A.  The  normal 
strain  at  point  A  is  a  function  of  L2i  only  and  is  given  by: 


,  -(L0.F-.+L1F_-HVL)r    F0 
e      (L  )     =     1f   2i  3  1  2  1  o  _2, 
eyAlL2i)     E1 I       +  AJ 


Therefore,  taking  the  difference  of  two  normal  strains  whose 
Y  coordinates  are  L2j_  =  L22  and  L2j_  =  L21,  such  that  L22  > 
L21'  yields 


(-L  9+L  ,  )r  F 


"v  *  22'    y^v  21 y  EI 

JA        ■'A 


Solving    for   F3    yields 


[cy    '42'-^    <L2l"EI 

F        =     _^ I* (1) 

3  (-L22+L21)ro 
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Figure  12.   Strain  Measurement  Locations 
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Alternatively,  using  the  same  rationale  for  point  B 
and  taking  the  difference  of  two  normal  strains  with 
coordinates  Y  =  L22  and  Y  =  L21»  yields 

„  ,     (L22-L2l'  ro  Fl 
SyB<L22>  "  eyB(L2l'  =  EI 

Solving  for  F]_  yields 

[fc  (L_J  -e.  (L01)]EI 
Yn     22    y   21 

F   =  5 E- (2) 

1         (L22  -L21}  ro 

At  points  A  and  C,  the  Mr  bending  moment  normal 
stress  components  are  equal  but  in  opposite  directions. 
Therefore  adding  the  normal  strains  at  these  points  yields 

2F2 
VL2l'  +  \     22       =  "^ 


Solving  for  F2  yields 


F2   =   -±* j-5 (3) 


With  the  three  load  forces  known,  the  moment  loads 
can  be  determined  directly  from  the  strain  equations. 
Solving  the  normal  strain  at  point  A  at  L2i  for  M^  yields 

F0   (L0,F0+L,F_)r 

1  2  o,  (4) 


Mj_  =  -[-EC   (L21)  +- j ] 


*A 
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Solving  the  normal  strain  at  point  B  at  L2i  for  M3  yields 


M3  =  J-[-E.   (La)  +4  +  Jiji^  (5) 

o    •'B 


Solving  the  shear  stress  at  point  B  at  L2i  for  M2  yields 


J  F3Q   LlFlr 

^  =  p*xM  '-k*^  (6) 

r    o 


7.   Previous  Load  Determinations 

As  stated  earlier,  the  previous  work  assumed  the 
principal  loads  were  F^  ,  F3 ,  and  M3  and  the  other  loads 
were  negligible.  The  further  assumption  that  F]_  did  not 
contribute  to  the  bending  moment  about  the  Z  axis  was  also 
made.   Thus  the  strain  equations  of  Table  I  reduce  to: 


<   (L  )   -   ~21  ' 3  ^ 

rA  -i]      EI 


-M  r 


1  "LlFlr   •    F3Q 
*y      (L21>      =     G[-f^  +  -k] 
ZB 


Solving   for   the    loads   yields 


F3     =    L^-eyJL2l)  <7> 

21  o     yA 


"3    =    =^8y  <L2!>  <8> 
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Vi  -  -  l^w  -  S1      '  '  <9> 

In  a  later  part  of  this  report,  the  results  of  equations 
(1)  ,  (5)  ,  (6)  and  equations  (7)  ,  (8)  ,  (9)  will  be  compared 
in  order  to  evaluate  their  predictions. 

C.   STRAIN  GAGE  PLACEMENT 

1.  Strain  Gage  Type  and  Location 

In  order  to  solve  equations  (1)  ,  (2)  ,  and  (3)  the 
axial  strain  at  five  locations  is  required.  Additionally, 
the  shear  strain  is  required  at  one  location  for  equation 
(6) .  Table  II  summarizes  the  strain  gage  types  and 
locations,  and  Figure  12  shows  their  locations. 

2 .  Determination  of  the  Shear  Stress 

Although  the  normal  strains  can  be  measured  directly 
from  the  gages,  additional  computations  are  required  to 
determine  the  shear  stress. 

Figure  13  depicts  a  45  degree  rectangular  rosette 
strain  gage.   The  applicable  strain  equation  is  [Ref.  3]: 


Y    =  2e   -  e  -  e 

yz      yz    z    y 


Substituting  this  expression  into  the  expression  for  M2 
yields 


v  o      L  F  r 


r 


yzg  21    zB  Zl         yB 
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TABLE  II 
STRAIN  GAGE  TYPES  AND  LOCATIONS 


Point 


Al 


Moment 
Arm 


L21 


Type 


Linear 


Quantity 
Measured 


A2 


Bl 


L22 
L21 


Linear 


4  5  degree     tz>*y>*yz 

rectangular 

rosette 


B2 


CI 


L22 
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Linear 


Linear 


e, 


Figure  13.   45  Degree  Rectangular  Rosette 
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Table  III  summarizes  the  force/moments  equations  for  this 
work  and  Table  IV  for  the  previous  study.  Appendix  B 
contains  the  Fortran  computer  code  for  these  equations. 

D.   CALIBRATION  PROCEDURES 
1.   Approach 

Static  loads  were  applied  to  the  support  structure 
to  evaluate  whether  the  analytical  expressions  of  Table  I 
correctly  predict  the  measured  strains  at  locations  A,  B, 
and  C.  Figures  14,  15,  16,  and  17  depict  the  methods  by 
which  the  static  loads  were  applied  to  the  instrumented 
drill  support  structure.  The  static  loads  consisted  of 
round  stock  cut  into  approximately  ten  pound  segments.  Each 
segment  weight  was  measured  using  an  Instron  Corporation 
Material  Testing  System  100  pound  test  cell.  For  this  test 
cell  the  accuracy  of  the  measurement  was  +•/-  1.0  %  of  the 
indicated  load,  or  0.02  %  of  the  test  cell  capacity  (100 
lbs)  . 

The  weight  holding  apparatus  depicted  in  the  static 
loading  figures  was  of  sufficient  strength  to  hold  the 
required  testing  loads.  However,  due  to  misalignment  in  the 
support  structure,  the  loads  were  not  applied  exactly  along 
the  desired  coordinate  directions   associated  with   the 
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TABLE    III 
LOAD   EQUATIONS 


tEy<L22'   "  £y <L21)!EI 
F,      =     • 


(L22  "L21)ro 


[£y   <L2l'   -£y   'L21)1EA 

F       = - 

r2  2 


IeyA(L22)   -£yA(L21)]EI 
F3     "  (-L22+L21)rQ 


_  F-        (L01F-.  +L..F  )r 

M       _       I  r  _        ,T      v  2            21  3       1  2     o1 

Ml     "     F" hE£ya(   2^    +  "A 1 3 

o         JA 


F3Q       L1F1r 


M,     =     -(G[2£  7    (L91)-£,    (L,,)-e      (L_)]    --^--+-±~-) 

z  —  yz        21        z„     21       v^     21  Jt  J 


'B  *  B 


.  F„        L01  F,  r 

o  ^B 


TABLE  IV 
LOAD  EQUATIONS  OF  PREVIOUS  WORK 


3     roL2i  ^ 

«3     =     "|M3 

Vl    =    -  I'^y^l'  "  %«4l>  "  £yB(L21)]   "  ¥' 
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Figure  14.   F2  Static  Loading 
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Figure  15.  F^   Static  Loading 
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Figure  16.   F3  Static  Loading 
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Figure  17.   Fx  and  F3  Static  Loading 
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support  structure.  Figure  18  illustrates  the  problem 
encountered  when  a  load  was  attempted  along  the  Z-axis.  To 
overcome  this  obstacle,  an  angle  position  indicator  was  used 
to  measure  the  angles  between  the  orthonormal  axes  and  the 
line  of  action  of  the  static  loads.  Only  two  of  the  three 
angles  were  measurable  and  the  following  relationship  was 
used  _to  determine  the  third  angle: 


1  =    \Jcos2et   +  cos2$  +  cos2S 


The  component  of  the  static  load  acting  on  each  of  the 
orthonormal  axes  was  then  determined  using  the  following 
relationships : 

F;l   =   F  cos  o( 

F2   =   F  cos  p 

F3   =   F  cos  S 

As  can  be  seen  from  Table  III,  each  of  the  forces  is 
equal  to  a  strain  quantity  times  a  constant,  i.e, 

Fl  "  t£yB(L22)  -  £yB(L21) ]*constant 

F2  =  C£yA(L2l)  +  £yc(L2i)  ]*constant 

F3   =    [^(^22-)    ~    eyA(L2i)]*constant 
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Figure  18.   Static  Force  Line  of  Action 


Therefore,   in   comparing   the   measured   versus   calculated 
values,  these  strain  quantities  were  used  to  determine  the 
accuracy  of  the  analytical  strain  equations. 
2 .   Static  Loading  Results 

Figure  19  shows  the  plot  of  [ev  (L22)~ey  (L2l)  3  * 
both  measured  and  calculated,  versus  an  applied  F^  static 
load  of  30.74  to  99.65  pounds.  As  can  be  seen  from  the 
figure,  there  is  a  negligible  difference  between  the 
measured  and  calculated  values.  Figures  20  and  21  show  the 
plots  of  the  F2  and  F3  component  forces  with  their 
respective  strain  quantities.  The  substantial  difference 
shown  in  the  later  figures  will  be  discussed  below. 
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Figure  22  shows  the  measured  versus  calculated 
strain  plot  for  an  F3  load  of  from  30.77  to  198.66  pounds. 
Again  the  measured  and  calculated  values  of  the  principal 
strain  quantity  correlate  well.  Figures  2  3  and  24,  for  the 
transverse  force  components,  exhibit  the  same  disporportion- 
ate  errors  as  was  found  with  the  F-^    static  load. 

Figures  25  and  2  6  show  the  plots  for  major  X  and  Z 
axis  loading.  This  loading  was  accomplished  by  placing  the 
support  fixture  at  an  angle  of  44.5  degrees  with  respect  to 
the  load  line  of  action.  Therefore,  a  static  load  of  30.77 
to  198.66  pounds  produced  F]^  and  F3  loads  of  22.83  to  147.42 
pounds  and  20.61  to  133.06  pounds  respectively.  Again  the 
strain  equations  produced  results  very  close  to  the  measured 
values  with  the  exception  of  the  transverse  F2  force,  as 
shown  in  Figure  27. 

The  results  of  applying  an  F2  static  load  are  shown 
in  Figure  28.  F2  loading  was  accomplished  by  placing 
weights  of  10.0  to  100.02  pounds  directly  on  top  of  the 
support  fixture.  In  this  test,  there  were  considerable 
differences  in  the  measured  and  calculated  values. 
3  .   Sources  of  Error 

The  errors  encountered  with  the  transverse  force 
components  and  F2  loading  are  attributed  to  the  nonsymmet- 
ric  shape  of  the  tube,  the  nonrigidity  of  the  tube,  and  the 
nonhomogeneous  tube  material.  These  factors  all  contribute 
to  change  the  centroidal  axis  location. 
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Figure  22.   F3  Load,  F3  Plot 
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Figure   25.      Fx   and   F3   Loads,    Fx   Plot 
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Figure  26.  T±   and  F3  Loads,  F3  Plot 
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Figure   28.      F2   Load,    F2   Plot 
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As  can  be  seen  in  Figures  6  and  7  the  location  of 
the  netitral  axis  for  the  bending  moments  coincides  with  the 
centroidal  axis.  Also,  the  magnitude  of  the  stress  is 
directly  proportional  to  the  distance  from  the  neutral  axis, 
and  the  stress  distribution  is  symmetrical  about  it.  These 
facts  operate  with  the  non-ideal,  real  tube  discussed  above 
to  create  a  changing  location  of  the  centroidal  axis  and, 
consequently  distort  the  stress  distribution.  As  shown  in 
the  transverse  force  component  plots,  the  change  in 
centroidal  axis  has  a  substantial  impact  on  the  stress 
distribution  and,  consequently,  the  reliability  of  the 
strain  equations  to  determine  the  loading.  But  due  to  the 
relative  magnitudes  shown  in  the  preceding  figures,  these 
errors  were  considered  negligible  in  predicting  the  dynamic 
loads  during  a  drilling  operation. 

Similarly,  the  stress  distribution  resulting  from  an 
F2  loading,  Figure  11,  was  nonuniform  with  a  change  in  the 
centroidal  axis  location,  yet  agreement  was  good. 
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IV.   DRILLING  OPERATION  AND  RESULTS 

A.   EXPERIMENTAL  PROCEDURES 

The  strain  gage  amplifier  assembly  and  strip  chart 
recorder  were  energized  and  given  a  one  hour  warm-up  period 
prior  to  commencement  of  recording  data.  The  air  compressor 
was  started  one-half-hour  prior  to  the  drilling  operations 
to  allow  sufficient  time  to  pressurize  the  accumlator  to  167 
psig.  After  the  one  hour  warm-up  period,  the  amplifiers 
were  calibrated  to  4.84  volts  per  10,000  micro-strains.  The 
penetrator  trigger  was  taped  in  the  full  open  position  and 
air  flow  was  controlled  by  an  in-line  valve.  The  penetrator 
was  positioned  via  the  milling  machine's  two  motors  to 
permit  the  penetrator  an  angle  of  attack  of  8.48  degrees  and 
a  three  diameter  separation  between  drill  hole  centers. 

At  the  commencement  of  a  drilling  operation,  the  in-line 
air  valve  was  fully  opened  and  the  Z-axis  motor  of  the 
milling  machine  energized  to  a  bed  speed  of  14  inches  per 
minute.  All  systems  were  stopped  after  the  penetrator 's 
drill  had  completed  drilling  the  hole. 

The  strip  .  chart  recordings  were  then  converted  to 
micro-strains.  A  sample  of  the  strip  chart  readings  is 
shown  in  Figure  29.  Since  the  conclusions  of  this  work  are 
to  be  used  in  the  design  of  the  penetrator  support 
structure,  a  worst-case  set  of  data  was   needed.   Therefore, 
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a  french  curve  was  used  to  curve-fit  the  maximum  values  of 
the  recorded  strains.  To  accurately  reflect  the  continous 
time  effects  of  the  drilling,  data  at  every  0.2  seconds  was 
used.  The  micro-strain  data  was  then  entered  into  the 
static  analysis  Fortran  program  of  Appendix  B,  which  is  the 
computer  code  equivalent  of  the  force/moment  equations  of 
Table  III  and  Table  IV. 

B.   AIR  SUPPLY  PRESSURE 

To  obtain  viable  data,  it  was  essential  that  the 
penetrator  successfully  complete  the  hole  drilling  without 
binding.  The  previous  work  used  an  operating  pressure  for 
the  penetrator  of  100  psig  [Ref.  2].  In  the  present  work, 
a  pressure  regulator  was  initially  installed  in  the  air  line 
to  reduce  the  accumulator  pressure  from  167  to  100  psig. 
With  100  psig  supply  pressure,  the  penetrator  was 
unsuccessful  in  completing  a  drilling  operation  without 
binding.  The  supply  pressure  was  increased  in  steps  of  25 
psig  from  100  to  150  psig  and  no  run  was  achieved  without 
experiencing  binding.  Therefore  the  pressure  regulator  was 
removed  from  the  system  and  runs  were  achieved  without 
binding  with  pressures  of  158  and  164  psig. 

As  drillings  were  accomplished,  the  accumulator  supply 
pressure  decreased  from  167  to  150  psig,  at  which  point  the 
compressor  would  recharge  the  accumalators.  During  the  ex- 
perimental process  it  was  noted  that  any  attempts  to  drill 
with  supply  pressures  of  less  than  158  psig  binding  occured. 
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C.   DRILLING  DATA  RESULTS 

Three  drilling  operations  were  conducted  with  supply- 
pressures  of  158,  164,  and  164  psig.  For  the  plots,  the 
drillings  are  differentiated  by  their  duration  times  as 
illustrated  in  Table  V. 

TABLE  V 
LOAD  RESULTS  DESIGNATION 


Run   Pressure   Duration 

No.     (PSIG)     Time  (sec)    Functions 

1  158  6.0       FifS.O) ,F2(6.0) ,F3(6.0) , 

Mi (6.0)  ,M2(6.0)  ,M3(6,.  0) 

2  164  7.2       F1(7.2)/F2(7.2) ,F3(7.2) 

M1(7.2),M2(7.2),M3(7.2) 

3  164  10.0       F1(10.0) ,F2 (10.0) ,F3(10.0) , 

M1(7.2) ,M2(7.2) ,M3(7.2) 


1 .   Forces 

The  force  loads  consist  of  a  thrust  component,  F3 , 
acting  along  the  axis  of  the  drill  and  two  transverse 
components,  F^  and  F2 ,  acting  in  a  plane  perpendicular  to 
the  axis  of  the  drill. 

As  can  be  seen  on  Figure  30,  the  force  F3  was  always 
in  the  positive  direction,  i.e.,  pushing  against  the  drill. 
The  graphs  indicate  the  values  of  F3  are  fairly  consistent 
no  matter  what  the  drilling  time  was.  For  instance,  the 
maximum  values  were  89.34,  88.47,  and  98.01  pounds  for  the 
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three  runs.  Also,  the  data  shows  that  F3  increases  to  a 
peak  as  the  drill  penetrates,  then  drops  in  magnitude  prior 
to  increasing  to  a  final  maximum  peak.  The  data  relates  to 
the  geometry  of  the  drill  tip  (Figure  31)  as  the  drill 
penetrates.  The  first  peak  is  caused  by  the  knuckle,  the 
second  by  the  trailing  edge. 

The  transverse  forces  on  the  other  hand,  vary  in 
direction  from  one  run  to  the  next  (Figures  32,  33).  The 
direction  change  is  attributed  to  the  movement  of  the  drill 
tip  prior  to  stablizing  in  a  drilling  configuration;  i.e., 
seating  into  the  test  panel.  If  the  initial  biting  movement 
of  the  drill  tip  is  referred  to  as  "walking,"  then  the 
direction  and  distance  of  the  walking  contribute  to  both  the 
magnitude  and  direction  of  the  transverse  loads.  Once  a 
stable  configuration  was  achieved  in  a  run  the  force 
direction  for  that  individual  run  did  not  change.  In  run  1, 
the  direction  of  the  drill  was  down  and  to  the  left  on  the 
surface  of  the  test  panel.  While  for  run  2,  the  direction 
was  up  and  to  the  left;  for  run  3,  the  direction  was  down 
and  to  the  right.  The  maximum  values  were  for  F^:  -24.28, 
-13.66,  and  44.67  pounds;  and  for  F2 :  -50.76,  42.35,  and 
-42.31  pounds.  The  maximum  values  for  F^  occured  at  the  end 
of  the  drill  run  and  for  F2  at  the  beginning  of  the  run. 
Therefore,  F2  was  the  governing  transverse  force  in 
determining  the  stablized  configuration.   The  F-±      force  was 
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Figure   30.      Drilling  Results   for  F3   Plot 
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Figure  32.   Drilling  Results  for  Fx  Plot 
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F2:  Y  AXIS  FORCE 
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Figure  33.   Drilling  Results  for  F2  Plot 
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characterized   by   increasing   magnitude   as   the   drill's 
trailing  edge  was  approached.   The  F2  force  exhibits  sharp, 
momentary  increases  in  magnitude  followed  by  relatively  long 
periods  at  constant  value. 
2.   Moments 

The  direction  and  magnitude  of  the  moments  were 
coupled  to  the  effects  of  the  transverse  load  seating 
conditions.  For  run  1,  the  M^  moment  about  the  X  axis  was 
positive  with  a  maximum  value  of  1137.88  inch-pounds  (Figure 
34),  and  for  run  2,  M^  was  negative  with  a  maximum  value  of 
-996.75  inch-pounds.  While  for  run  3,  the  value  of  M-^ 
varied  from  -37.53  to  1003.91  inch  pounds.  Each  time,  the 
maximum  values  occured  at  the  beginning  of  the  drilling. 

The  M2  moment  about  the  Y  axis  oscillated  between 
negative  and  positive  values  in  all  three  runs  (Figure  35) . 
The  ranges  of  values  were  as  follows: 

Run  1:   -427.4  5  to  72.2  8  inch-pounds 
Run  2:   -147.11  to  47.99  inch-pounds 
Run  3:   -309.63  to  803.71  inch-pounds 
The  peak  values  were  at  the  trailing  edge  of  the  drill  bit. 
The   M3   moments   about   the   Z   axis   were   also 
oscillatory,  but  were,  for  the  most  part,  positive  in  nature 
(Figure  36)  .   This  positive  nature  correlates  to  the  fact 
that  binding  did  not  occur  for  these  runs.   The  range  of 
values  obtained  were: 
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Ml:  X  AXIS  MOMENT 
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Figure  34.   Drilling  Results  for  Mx  Plot 
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M2:  Y  AXIS  MOMENT 
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Figure  35.   Drilling  Results  for  M2  Plot 
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M3:  Z  AXIS  MOMENT 
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Figure  36.   Drilling  Results  for  M3  Plot 
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Run  1:   -3.37  to  63.23  inch-pounds 
Run  2:    0.0  to  67.55  inch-pounds 
Run  3:   -18.03  to  160.52  inch-pounds 

As  with  M2/  the  peak  values  were  at  the  drill's  trailing 

edge. 

D.  COMPARISON  WITH  PREVIOUS  WORK 

Figures  37  through  4  5  are  plots  which  compare  the 
results  of  this  work  to  the  forces/moments  computed  using 
the  equations  of  the  previous  work  in  Table  IV.  In  the 
plots,  various  quantities  are  differentiated  by  a  PW 
preceeding  values  of  the  previous  work.  The  values  of  F3 , 
L1*F1'  and  M3  correlated  fairly  well  for  the  two  studies  and 
thus  complement  each  other's  validity. 

The  differences  in  the  F3  force  are  due  mainly  to  the 
assumption  in  the  previous  work  omitting  the  M^  and  L^*F2 
bending  moments  which  contribute  to  the  state  of  strain  at 
point  A. 

Similarly,  the  twisting  moment  at  point  B  is  not  only  a 
function  L1*F1  but  also  the  moment  M2  which,  depending  on 
its  direction,  adds  or  subtracts  from  L^*F1.  The  normal 
strain  at  point  B,  is  not  only  a  function  of  M3 ,  but  also 
the  bending  moment  caused  by  L2*F]_. 

E.  ERROR  ANALYSIS 

For  the  strip  charts,  a  one-half  of  a  division  error 
equates  to  the  following  accuracies  for  the  forces: 
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F3:  Z  AXIS  FORCE  (M1=M2=F2=0) 
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Figure   37.      Current  Versus   Previous  Work:      F3    Plot   1 
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F3:  Z  AXIS  FORCE  (M1=M2=F2=0) 
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Figure   38.      Current  Versus   Previous  Work:      F3   Plot   2 
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F3:  Z  AXIS  FORCE  (M1=M2=F2=0) 
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Figure  39.   Current  Versus  Previous  Work:   F->  Plot  3 
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L1*F1:Y  AXIS  M0MENT(M1=M2=F2=0) 
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Figure  40.   Current  Versus  Previous  Work:   Mn  Plot  1 
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Ll*Fl:Y  AXIS  MOMENT(M1=M2=F2=0) 
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Figure  41.   Current  Versus  Previous  Work:   M2  Plot  2 
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L1*F1:Y  AXIS  M0MENT(M1=M2=F2=0) 
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Figure  42.   Current  Versus  Previous  Work:   M2  Plot  3 
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M3:  Z  AXIS  MOMENT(M1=M2=F2=0) 
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Figure  44.   Current  Versus  Previous  Work:   M3  Plot  2 
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M3:  Z  AXIS  MOMENT(M1=M2=F2=0) 
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Figure   43.      Current  Versus   Previous  Work:      M3    Plot   1 
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M3:  Z  AXIS  MOMENT(M1=M2=F2=0) 
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Figure   45.      Current  Versus   Previous  Work:      M3   Plot   3 
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Fl:  +/-  2.17  pounds 

F2:  +/~  42.31  pounds 

F3:  +/-  8.67  pounds 
The  results  of  the  static  tests  of  the  equations  yielded 

the  following  maximum  error  percentages  for  measured  versus 
calculated  values: 

Fl:  +/-   0.4  % 

F2:  +/"  10.2  % 

F3:  +/"  1.1  % 
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V.   CONCLUSIONS 

An  investigation  of  the  requirements  for  a  robot  arm  to 
manipulate  an  aircraft  skin  penetrator/applicator  was 
conducted.  The  penetrator  was  a  USAF-developed  device  which 
was  designed  to  penetrate  the  fuselage  of  an  aircraft  and 
inject  a  gaseous  fire  fighting  agent  (Halon)  into  the 
interior  of  the  aircraft  to  extinguish  fires.  Experiments 
were  designed  to  identify  worst-case  drilling  loads,  and  a 
test  apparatus  was  used  to  mimic  the  drilling.  Test  data 
were  collected  on  the  effects  of  drilling  through  sheet 
aluminum  alloy.  Forces  and  moments  were  measured  by  strain 
gages  during  actual  drilling. 

The  test  data  indicates  that  the  loads  at  the  drill  tip 
are  significantly  influenced  by  drill  walking  prior  to 
stabilizing  and  penetrating  of  the  drilling  surface.  This 
is  attributed  to  significant  transverse  load  forces  (F^  and 
F2)  which  act  in  a  plane  perpendicular  to  the  drill  axis. 
These  transverse  forces  act  with  the  length  of  the 
penetrator  to  exert  significant  moments  on  the  support 
structure.  Significant  factors  in  the  walking  phenomenon 
include:  the  age  and  condition  of  the  drill  tip;  the  type 
of  material  to  be  drilled;  and  the  angle  of  attack  between 
the  drill  axis  and  the  drilling  surface,  all  of  which 
contribute  to  the  direction  and  the  magnitude  of  the  drill 


156 


walking  prior  to  penetration.  The  drill  walking  in  turn 
determines  the  directions  and  magnitudes  of  the  transient 
loads. 

The  following  is  a  summary  of  the  maximum  forces  and 
moments  for  each  of  the  runs  recorded: 

Run  1  Run  2       Run  3 

F1    (lbs)        -24.28  -13.66        44.67 

F2  (lbs)        -50.76  42.35       -42.31 

F3  (lbs)         89.34  88.47        98.01 

Mx  (in-lbs)  1137.88  -996.75  1003.91 
M2  (in-lbs)  -427.45  -147.11  803.71 
M3  (in-lbs)         63.23  67.55       160.52 

From  this  summary,  a  negative  F2  results  in  a  positive  K-^ 
moment  on  the  drill  tip  where  the  reverse  case  is  true  for 
positive  F2 .  As  for  F^,  its  direction  effects  the  resulting 
moment  M2 .  The  direction  of  M2  is  positive  when  F^  is 
positive  and  vice  versa.  With  regard  to  M3 ,  F]_  is  the 
dominant  transverse  force,  but  not  to  the  extent  it  was  for 
M2,  this  is  true  because  the  moment  arm  in  the  M3  case  is 
the  distance  from  the  centerline  of  the  drill  to  a  point  on 
the  support  structure  which  is  much  less  than  the  length  of 
the  penetrator  in  the  M2  case.  But  as  can  be  seen  above, 
the  magnitude  of  M3  is  significantly  higher  for  run  3  where 
a  positive  F^  was  realized  vice  the  other  runs  when  a 
negative  F^  resulted  from  the  drill  walking. 
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With  regard  to  the  previous  study  of  Reference  2,  the 
effects  of  drill  wear  are  an  increased  axial  force  and/or 
torque  requirement.  The  location  and  nature  of  the  wear 
governs  which  symptom  dominates.  Chisel  edge  wear  results 
in  increased  axial  force,  land  wear  results  in  increased 
torque,  and  wear  at  the  corners  results  in  increases  of  both 
[Ref.  4].  The  thrust  (F3)  and  torque  (M3)  of  the  previous 
work  for  an  8  degree  angle  of  attack  were  of  the  magnitudes, 
respectively,  of  280  pounds  and  475  inch-pounds  for  100  psig 
supply  air,  and  205  pounds  and  255  inch-pounds  for  140  psig 
supply.  Therefore,  the  combination  of  the  different  wear 
effects  discussed  above  combined  with  an  increase  of  supply 
air  pressure  to  160  psig,  would  explain  the  higher  thrust 
and  torque  required  for  the  worn  drill  used  in  the  previous 
study  and  the  lower  requirements  determined  in  this  study. 
This  was  also  visually  evident  in  the  resulting  sharp,  clean 
holes  for  the  current  study  and  the  rough,  jagged  holes 
drilled  with  the  worn  drill.  Another  contributing  factor  to 
the  differences  in  the  torque  moment,  is  the  fact  that  the 
previous  study  did  not  deduct  the  bending  effects  caused  by 
Fi  times  the  length  from  the  centerline  of  the  drill  to  the 
strain  gage  location.  Similarly,  the  previous  assumption 
that  there  was  no  M2  at  the  drill  tip,  just  twisting  of  the 
support  structure  caused  by  F]^  times  the  length  of  the 
penetrator,  disallowed  any  correlation  between  the  two 
studies  with  respect  to  F^  and  M2 . 
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VI.   RECOMMENDATIONS  FOR  FUTURE  WORK 

Further  investigation  for  two  specific  areas  are 
required.  The  first  is  the  phenomenon  of  the  drill  walking 
and  the  resulting  effects  it  has  on  the  drill  stabilizing 
and  penetrating.  The  drill  has  a  possible  choice  of  four 
quadrants  in  which  to  move.  These  would  result  from  the 
following  combinations  of  transverse  forces: 

Fi   and  F2 

Fi   and  -F2 

-F^   and  F2 

-F]_   and  -F2 

Drilling  data  should  be  obtained  for  all  the  combinations 

and  analyzed  for  drill  support  effects  and  repeatability. 

The  second  area  of  concern  is  the  life  cycle  of  the 
penetrator  drill.  The  thrust  and  torque  requirements  are 
both  a  function  of  drill  wear.  In  view  of  the  fire  fighting 
mission  of  the  device,  a  definitive  useful  life  cycle  for 
the  drill  bit  should  be  identified.   -  -  . 

With  regard  to  the  experimental  procedures,  the 
methodology  of  this  study  is  sufficient  but  improved  data 
recording  is  required  for  further  work.  A  process  for 
directly  converting  the  continuous  time  output  of  the  strain 
gages  into  an  input  for  the  strain-to-load  conversion 
computer  program  should  be  developed. 
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APPENDIX  A 


STRAIN  EQUATIONS  FORTRAN  COMPUTER  CODE 

x  LAWRENCE,  D.A. 

DRILL1  PROGRAM  CLACULATES  THE  RESULTING  STRAINS  GIVEN  THE  APPLIED 
MOMENTS  AND  FORCES  AT  THE  DRILL  TIP. 
VARIABLES  USED  ARE: 


TH 
RO 
RI 

E: 


RADIAL  THICKNESS  (IN) 
OUTSIDE  CYLINDER  DIAMETER  (IN) 
INSIDE  CYLINDER  DIAMETER  (IN) 
YOUNGS  MODULUS  OF  ELASTICITY  (LBS/IN**2) 


POISON  RATIO 


POI 

PI 

G:  SHEAR  MODULUS 

D:  MEAN  DIAMETER 


OF  ELASTICITY 
(IN) 


(LBS/INXX2) 


L21 
L22 
LI: 
Fli 
F2. 

F3: 
Ml  : 
M2: 
M3: 

EYA21 
EYA22 
EYB21 
EYB22 
EYC21 

:yzbzi 

EZB21: 
EYZB21 


DISTANCE 
DISTANCE 


FROM 
FROM 


CENTERLINE 
CENTERLINE 


OF 
OF 


DRILL 
DRILL 


TO 
TO 


STRAIN 
STRAIN 


DISTANCE  FROM  DRILL 
X  AXIS  LOAD  (LBS) 

LOAD 

LOAD 

LOAD 

LOAD 

LOAD 


TIP  TO  CYLINDER  CENTROID 


GAGE  (IN) 

GAGE(IN) 

(IN) 


Y  AXIS 

Z  AXIS 

MOMENT 

MOMENT 

MOMENT 
NORMAL 
NORMAL 
NORMAL 
NORMAL 
NORMAL 
SHEAR 
NORAML 


(LBS) 
(LBS) 
ABOUT 
ABOUT 
ABOUT 
STRAIN 
STRAIN 
STRAIN 
STRAIN 
STRAIN 
STRAIN 
STRAIN 


POINT 
POINT 
POINT 
POINT 
POINT 
POINT 
Y-Z  PLANE  STRAIN  ON  45 


X 

Y 

Z 
AT 
AT 
AT 
AT 
AT 
AT 
AT 


AXIS  (IN-LBS) 
AXIS  (IN-LBS) 
AXIS  (IN-LBS) 
POINT  A  DISTANCE 


DISTANCE 
DISTANCE 
DISTANCE 
DISTANCE 
DISTANCE 
DISTANCE 
DEGREE  RO: 


L21 

L22 
L21 
L22 
L21 
L21 
L21 
ETTE 


LEG  DISTANCE  L21 


X( I) CONSTANTS  INCLUDE  I , Q, RO, D. P I , G, E 
A(I.J):  REPEATED  CALCULATIONS  IN  EQUATIONS 


*XXX»X<i*XX)t<XX**»XX*<X<*KX*X*»«»XX»»>i(XX.<XX*»i(J(*XXX»»i(X»»J()<XXi()<X»»XXX) 


DOUBLE 


PRECISION  RO 
A22,B21,B22. 
F1,F2,F3,M1, 

J 
PARAMETERS 


RI,TH,E,POI,PI,G,D,L21,122,Ll 
EYA21,EYA22,EYB21,EYB22,EYC21 
M2,M3,Y1,Y2,Y3 


+  A1.A21 
+  SYZ321 

INTEGER 
*DEFINE  KNOWN 

R0=0.5 

TH=0.1 

RI=RO-TH 

£  =  29000000. 0 

POI=0.29 

PI=4.0*ATAN(1.0) 

G  =  E/(2.0*(1  .0  +  POD) 

D=(2.0*RO)-TH 

L21=2.424 

L22=6 .424 

L1=0.0 

^CALCULATE  CONSTANT  TERMS 

X1=PI*THXD*E 
X2=4.0XRO/(PIXE*(ROXX4-RIXX4)) 

X3=8.0X(R0XX3-Rlxx3)/(3.0XGXPIX(THXX2)X(DXX3)) 
X4=GXPIXTH*(DXX2) 
XINPUT  THE  LOADS 
DO  10  J=l,4 
READ(14,x)  F1,F2,F3,M1,M2,M3 


X1.X2,X3.X4 
EYZB21.E2B21 


IN  EQUATIONS 
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A1=F2/X1 

A21=-(L21*F3+L1XF2+M1)XX2 
A22=-(L22XF3+L1XF2+M1)XX2 
B21=(L21*F1-M3)XX2 
B22=(L22xFl-M3)*X2 
XCALCULATE  THE  STRAINS 
EYA21=(A21+A1) 
EYA22=(A22+A1) 
EYB21=(B21+A1) 
EYB22=(B22+A1) 
EYC21=(-A21+A1) 

SYZB21=(-2.0X(L1XF1-M2)/X^)+(X3XF3) 
£ZB21=-P0I*EYB21 
EYZB21=(EYB21+EZB21+SYZB21)/2.0 
XPRINT  LOADS  FOR  INPUT  TO  EASYPLOT 

0PEN(UNIT=21,FILE=,DATA1,,STATUS=,0LD») 

WRITE(21,5) 
5     F0RMAT(T5,,F1ST13,,F2,,T21,,F3VT29,,M1,/T37,'M2ST^5,,M3V) 

WRITEC21,15)  F1,F2,F3,M1,M2,M3 
XPRINT  OUTPUTS  OF  STRAINS  AND  LOADS 
15     FORMATC6F10.2//) 

WRITEC21,25) 
25     F0RMATCT5,  '  EYA21 » ,T18,  •EYA22S/) 

WRITE(21,35)EYA21,EYA22 
35     F0RMAT(2E12.<+/) 

HRITEC21,45) 
45     F0RMATCT5,  •  EYB21  •  ,T18 ,  '  EYB22*  ,T30,  » EYZB21  ■ ,  T42,  f  EZB21 '/ ) 

HRITE(21,55)  EYB21 , EYB22, EYZB21, EZB21 
55     FORMATS  El  2.4/) 

WRITE<21,65) 
65     F0RMATCT5,  'EYC21V) 

HRITE(21,75)  EYC21 
75     F0RMATCE12.4//72C »+•)//) 
*PRINT  STRAINS  TO  FILE  FOR  CALIBRATION  PCOTS 

EYA21=1000000*EYA21 

EYA22=1000000XEYA22 

EYB21=1000000*EYB21 

EYB22=1000000*EYB22 

EYZB21=1000000*EYZB21 

EZB21=1000000*EZB21 

EYC21=1000000*EYC21 

WRITEC20 , X)  EYA21 , EYA22, EYB21 , EYB22, EYZB21 , EZB21 , EYC21 

Y1=EYB22-EYB21 

Y2=EYA21+EYC21 

Y3=EYA22-EYA21 

0PEN(UNIT=22,FILE='DATA2,,STATUS=,0LD») 

0PEN(UNIT=24,FILE=,DATA3,,STATUS=,0LD,) 

0PEN(UNIT=25,FILE=,DATA4f,STATUS=,0LD») 

HRITEC22,85)  F1,Y1 

WRITE(24,85)  F2,Y2 

WRITE(25,85)F3,Y3 
85     FORMAT(2F10.2) 
10    CONTINUE 

STOP 

END 
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APPENDIX  B 
LOAD  EQUATIONS  FORTRAN  COMPUTER  CODE 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXJ 

x  LAWRENCE,  D.A. 

THIS  PROGRAM  CALCULATES  THE  THE  LOADS  AT  THE  PENETRATOR  DRILL  TIP 
BASED  ON  MICRO-STRAIN  INPUTS.  ALSO  THE  PROGRAM  CALCULATES  THE  LOADS 
USING  THE  ASSUMPTIONS  THAT  F2,M1,AND  M2  ARE  ZERO  MADE  IN  THE  PREVIOUS 
WORK. 


VARIABLES  USED  ARE: 
T:TIME  (SEC) 

TH:  RADIAL  THICKNESS  (IN) 
RO:  OUTSIDE  CYLINDER  DIAMETER  (IN) 
I j  MOMENT  OF  INERTIA  (IN**4) 
E:  YOUNGS  MODULUS  OF  ELASTICITY 
POI:  POISON  RATIO 
PI 

G«  SHEAR  MODULUS  OF  ELASTICITY  CLBS/INXX2) 
D:  MEAN  DIAMETER  (IN) 

L21:DISTANCE  FROM  CENTERLINE  OF  DRILL 
L22:DISTANCE  FROM  CENTERLINE  OF  DRILL 


(LBS/INXX2) 


TO 

TO 


STRAIN 
STRAIN 


LI: DISTANCE  FROM  DRILL  TIP 
Fl:  X  AXIS  LOAD  FORCE(LBS) 
F2:  Y  AXIS  LOAD  FORCE(LBS) 
F3:  Z  AXIS  LOAD  FORCE(LBS) 
Ml:  MOMENT  LOAD  ABOUT  X 
M2«  MOMENT  LOAD  ABOUT  Y 
M3:  MOMENT  LOAD  ABOUT  Z 
EYA21 :  NORMAL  STRAIN  AT 
EYA22:  NORMAL  STRAIN  AT  POINT 
EYB21:  NORMAL  STRAIN  AT  POINT 
EY322:  NORMAL  STRAIN  AT  POINT 
CYC21:  NORMAL  STRAIN  AT  POWIT 
SYZB21;  SHEAR  STRAIN  AT  POINT 
EZB21:  NORMAL  STRAIN  AT  POINT 
EYZB21:Y-Z  PLANE  STRAIN  ON  45 
X(I):  CONSTANTS  INCLUDE  I,Q,RO 


TO  CYLINDER  CENTROID 


GAGE 

GAGE 
(IN) 


(IN) 
(IN) 


AXIS(IN-LBS) 
AXIS(IN-LBS) 
AXIS(IN-LBS) 
POINT  A  DISTANCE 


DISTANCE 
DISTANCE 
DISTANCE 
DISTANCE 
B  DISTANCE 
B  DISTANCE 
DEGREE  LEG 
,D,PI,G,E 


L21 
L22 

L21 
L22 
L21 
L21 
L21 
ROSETTE 


DISTANCE  L21 


A(I,  J)  REPEATED  CALCULATIONS  IN  EQUATIONS 
YF3:F3  LOAD  EQUATION  OF  PREVIOUS  WORK  (LBS) 
YM3:M3  LOAD  EQUATION  OF  PREVIOUS  WORK  (LBS) 
YM2:M2  LOAD  EQUATION  OF  PREVIOUS  WORK  (LBS) 

XXXXXXXXXXXXXXXXXXXXX***X*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DOUBLE  PRECISION  RO . RI , TH , E, POI , PI , G, D. L21 , L22, LI , XI , X2 , X3 , 
+  A1,A21,A22,B21,B22,EYA21,EYA22,EYB21,EYB22,EYC21,SYZB21,EZB21, 
+  EYZB21,F1,F2,F3,M1,M2,M3,T 

INTEGER  J 
XDEFINE  KNOWN  PARAMETERS 

R0=0.5 

TH=0.1 

RI=RO-TH 

E=29000000.0 

POI=0.29 

PI=4.0XATAN(1 .0) 

G  =  E/(2.0X(1  .0  +  POD) 

D=(2.0XRO)-TH 

L21=2.75 

L22=6.75 

Ll=23.125 
^CALCULATE  CONSTANT  TERMS 

X1=PIXTHXDXE 


IN  EQUATIONS 
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X2=4 . OXRO/(PIXEX(ROXX4-Rlxx«) ) 

X3=8.0X(ROxx3-Rlxx3)/(3.0xGxpix(THXX2)x(Dxx3)) 

X4=GXPIXTHX(DXX2) 

T  =  0.0 

DO  10  J=l,51 
XINPUT  MICRO-STRAIN  VALUES 

READ(8,X)  EYA21,EYA22,EYB21,EYB22,EYZB21,EZB21,EYC21 

EYA21=.000001XEYA21 

EYA22=.000001XEYA22 

EYB21=.000001XEYB21 

EYB22=.00000lxEYB22 

EYZB21=.000001XEYZB21 

EZB21=.000001XEZB21 

EYC21=.000001XEYC21 
XCALCULATE  FORCE  AND  MOMENT  LOADS 

F1  =  (EYB22-EYB21)/(X2X(L22-L2D) 
F2=Xlx(EYA21+EYC21)/2.0 
F3  =  (EYA22-EYA21)/(X2X(-L22+L2D) 
A1=F2/X1 

A21=(L21XF3+L1XF2)XX2 
A22=(L22XF3+L1XF2)XX2 
B21=-(L21XF1)XX2 
B22=-(L22XF1)XX2 
M1=(CA1-EYA21)/X2)-L21XF3-L1XF2 
M3=((A1-EYB21)/X2)+L21XF1 
SYZB21=2.0xEYZB21-EZB21-EYB21 
M2=(X4X(SYZB21-X3xF3)/2.0)+LlxFl 
XPRINT  LOADS  TO  FILE  FOR  INPUT  TO  EASYPLOT 
WRITE(14,2)  T,F1,F2,F3,M1,M2,M3 

2  FORMATC7F10.2) 

XCALCULATE  LOADS  BASED  ON  PREVIOUS  WORK  ASSUMPTIONS 

YF3=-EYA21/(X2XL21) 

YM3=-EYB21/X2 

YM2=X4x(X3*F3-SYZB21)/2 
XPRINT  OLD  AND  NEW  LOADS  OUANITIES  TO  FILE  FOR  EASYPLOT 

WRITE(15,3)  T,F3,YF3,M2,YM2,M3,YM3 

3  FORMATC7F10.2) 

xpRINTED  OUTPUT  OF  STRAINS  AND  LOADS 

0PErl(UNIT=23,FILE=lDATA5,,STATUS=,0LD,) 

WRITEC23.5)" 
5      F0RMATCT5, ,EYA21,,T1S, 'EYA22'/) 

WRITE<23,15)  EYA21,EYA22 
15     F0RMAT(2E12.4/) 

HRITE(23,25) 
25     F0RMATCT5, ,EYB21',T18, 'EYB22',T30, 'EYZB21 '  ,T42, 'EZBZIV) 

WRITE(23,35)  EYB21 , EYB22, EYZB21 , EZB21 
35     F0RMAT(4E12.4/) 

"  •  WRITE(23,45) 
45     F0RMAT(T5, 'EYC21'/) 

WRITE(23,55)  EYC21 
55     F0RMATCE12.4/) 

WRITE(23,65) 
65     F0RMAT('T»,T5,,F1,,T13,,F2,,T21,,F3,,T29,'M1,,T37,,M2,,T45,,M3V 
+  48CK*)/) 

WRITEC23,75)  T, Fl, F2, F3,M1,M2,M3 
75     FORMAT(7F10.2/72(»+')//) 

t=T+0.2 
10     CONTINUE 

STOP 

END 
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